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TECHNICAL  SUMMARY 

The  project’s  primary  objective  was  to  develop  suitable  theoretical  methods  to  simulate  the  transport  properties  of 
individual  molecules  and  provide  information  about  the  mechanisms  that  control  these  properties  in  the  context  of  efforts 
to  develop  “moletronic”  devices  and  systems  by  other  teams  in  the  Moletronics  Program.  The  task  was  accomplished: 
state-of-the-art  methods  were  developed  that  allow  the  calculation  of  the  current  voltage  characteristics  of  individual 
molecules  and  include  current-induced  atomic  relaxations  and  motions;  the  methods  were  employed  for  a  series  of 
molecules  that  were  fabricated  and  measured  by  the  team  headed  by  Mark  Read  and  James  Tour.  The  experimental 
observations  were  accounted  for  and  the  underlying  mechanisms  that  control  the  properties  of  the  molecule  were 
identified.  In  addition,  three-terminal  devices,  not  yet  fabricated  but  necessary  for  functional  systems,  were  simulated 
and  predictions  were  made.  The  results  were  published  in  five  articles  in  premier  journals  and  several  Conference 
Proceedings.  A  summary  of  the  results  is  given  in  this  report  and  the  original  articles  are  attached  as  Appendices. 

INTRODUCTION 

Silicon-based  microelectronics  is  reaching  the  level  of  miniaturization  where  quantum  phenomena  such  as 
tunneling  cannot  be  avoided  and  the  control  of  doping  in  ultrasmall  regions  becomes  problematical.  Though  it  is  likely 
that  silicon-based  technology  will  simply  move  to  a  different  paradigm  and  continue  taking  advantage  of  the  existing 
vast  infrastructure  and  manufacturing  capabilities,  novel  and  alternative  approaches  may  give  new  insights  and 
ultimately  may  usher  a  new  era  in  nanoelectronics.  Molecules  as  individual  active  devices  are  obvious  candidates  for  the 
ultimate  ultrasmall  components  in  nanoelectronics.  Though  the  idea  has  been  around  for  more  than  two  decades,1  only 
recently  measurements  of  current- voltage  characteristics  of  individual  molecules  have  been  feasible.  There  have  been 
significant  breakthroughs  in  experimental  measurements  that  offer  considerable  promise.2'7 

Methods  for  the  calculation  of  current  in  small  structures  placed  between  two  metal  electrodes  have  been  developed 
over  the  years,  but  actual  implementations  have  been  scarce.  For  molecules,  semiempirical  methods  have  been  used  to 


study  the  dependence  of  current  on  various  aspects  of  the  problem,8  but  quantitative  predictions  for  direct  comparison 
with  data  are  not  possible  because  values  of  parameters  under  current  conditions  cannot  be  determined  independently. 

In  the  1980's,  Lang  developed  a  practical  method  to  calculate  transport  in  the  context  of  imaging  atoms  with 
scanning  tunneling  microscopy.9  The  method  has  all  the  ingredients  needed  to  compute  current-voltage  (I-V) 
characteristics  of  single  molecules.  We  adopted  this  method  and,  in  addition,  we  developed  a  suitable  Hellmann- 
Feynman  theorem  for  the  calculation  of  current-induced  forces  on  atoms  that  allows  us  to  study  the  effect  of  current  on 
relaxations  and  ultimately  the  breakdown  of  molecules.10  With  these  tools,  we  have  carried  out  extensive  studies  of 
transport  in  molecules  whose  core  is  a  single  benzene  ring.  Such  molecules  have  been  synthesized  and  measured  by 
Reed  and  coworkers.2,3  Here,  we  summarize  the  most  important  results  of  the  recent  work.  The  method  of  calculation 
and  more  details  of  the  results  can  be  found  in  the  original  papers.1 1-14 

TRANSPORT  IN  A  SINGLE  BENZENE  RING 

Most  of  our  theoretical  work  has  been  carried  out  on  single  benzene  rings  because  the  small  number  of  atoms 
makes  the  calculations  practical.  Experimental  data  are  available  for  two-terminal  configurations. 

Fig.  1  shows  schematics  of  a  single  benzene  ring  bridging  the  gap  between  macroscopic  gold  electrodes.  A  sulfur 
atom  at  each  end  joins  the  benzene  ring  to  the  electrodes.  The  experimental  I-V  characteristic  is  shown  in  the  top  panel. 
The  middle  panel  shows  the  theoretical  results.11  We  will  discuss  the  third  panel  shortly. 

It  is  clear  from  the  figure  that  the  theory  reproduces  the  shape  of  the  I-V  curve  quite  well,  but  the  absolute 
magnitude  of  the  current  is  off  by  more  than  two  orders  of  magnitude.  We  address  each  of  these  issues  separately. 

In  Fig.  2  we  show  the  density  of  states  of  the  molecule  for  three  different  voltages:  0.01  V,  2.4  V  and  4.4  V  and 
mark  out  the  energy  window  between  the  left  and  right  quasi  Fermi  levels.  States  within  this  window  contribute  to 
transport.  We  see  that  there  is  virtually  no  density  of  states  in  the  small  window  at  small  voltages,  in  agreement  with  the 
slow  initial  rise  of  the  I-V  curve.  At  2.4  V,  the  ti*  states  of  the  molecule  enter  the  transport  window  and  give  rise  to  the 
first  peak  in  the  spectrum.  At  4.4  V,  the  n  states  of  the  molecule  enter  the  window  while  the  71*  continue  to  participate, 
giving  rise  to  the  second  peak  in  the  spectrum.  The  peak  at  2.4  V  is  somewhat  more  pronounced  in  the  theoretical  curve. 
The  observed  smoothing  is  likely  to  be  caused  by  interactions  between  the  electrons  and  vibrational  modes. 

In  order  to  explore  the  mechanism  that  controls  the  absolute  magnitude  of  the  current  we  performed  calculations  by 
inserting  an  extra  gold  atom  between  the  sulfur  atom  and  the  macroscopic  electrode  at  each  end  of  the  molecule  as 


shown  in  the  lower  part  of  Fig.  1.  There  was  a  dramatic  decrease  in  the  current,  bringing  its  value  much  closer  to  the 
experimental  value  (lower  panel  in  Fig.  1).  The  decrease  is  attributed  to  the  fact  that  gold  atoms  have  only  one  5  electron 
available  for  transport  and  s  electrons  do  not  couple  with  the  n  electrons  of  the  molecule.  The  gold  atoms  act  as  a 
quantum  mechanical  constriction.  To  test  the  idea  we  performed  calculations  by  replacing  the  gold  atoms  with  aluminum 
atoms.  The  latter  have  p  electrons  that  should  couple  well  with  the  n  electrons  of  the  molecule.  Indeed,  the  current 
jumped  to  its  initial  value.  An  additional  test  was  carried  out  with  three  gold  atoms  instead  of  a  single  gold  atom.  The 
current  was  again  at  its  full  value  because  the  three  s  orbitals  on  the  three  gold  atoms  can  form  the  appropriate  linear 
combinations  to  produce  sufficient  coupling. 

It  is  clear  from  the  above  that  molecules  determine  the  shape  of  the  I-V  characteristic,  but  the  nature  of  individual 
atoms  at  the  molecule-electrode  contact  determines  the  absolute  magnitude  of  the  current.  The  results  illustrate  the  power 
of  theory  to  contribute  to  device  design  ,  especially  “contact  engineering”. 

A  THREE-TERMINAL  DEVICE 

Transistor-like  behavior  in  molecular  three-terminal  devices  (molecular  transistors)  has  been  recently  demonstrated 
in  a  single  C60  bucky  ball  connected  to  gold  electrodes,4  in  organic  self-assembled  monolayers,5  and  in  single  and  multi- 
walled  carbon  nanotubes.6,7  The  source-drain  current  in  these  systems  (three-terminal  devices)  is  controlled  by  a  gate 
voltage. 

As  an  example  of  a  three-terminal  device  we  considered  the  same  benzene  ring  as  above  but  put  a  third  “terminal” 
in  the  form  of  a  capacitive  gate:  an  external  electric  field  across  the  molecule  produces  polarization  that  affects  the 
current  through  the  molecule.  Figure  3  shows  that  at  a  very  small  (0.01  V)  source-drain  voltage,  the  current  is  a  strong 
function  of  the  gate  voltage.  In  particular,  we  demonstrate  substantial  gain  at  the  resonant  peak. 

Recently,  we  explored  the  effect  of  symmetry  on  the  current-gate-field  (I-E)  characteristics  of  model  molecular 
transistors.  In  particular,  we  investigated  the  role  of  intrinsic  dipole  moments  on  the  current  modulation  in  these 
structures.  Figure  3  shows  the  two  core  molecules  that  have  been  constructed  from  the  benzene- 1 ,4-dithiolate  molecule 
by  substituting  one  of  two  hydrogen  molecules  by  hydroxile  groups  (OH).  The  transport  properties  in  these  two 
molecules  were  found  to  be  very  different.  The  amplification  of  the  current  at  the  resonant  tunneling  regime  is  a  few 
times  larger  in  the  asymmetric  molecule  which  is  a  subject  to  a  linear  Stark  effect  in  the  gate  electric  field.  At  the  same 
time,  resonant  tunneling  occurs  at  much  lower  voltage  than  in  the  symmetric  molecule. 


The  shape  of  the  I-E  characteristics  depends  on  the  value  of  the  source-drain  voltage  in  a  nonlinear  way,  thus 
manifesting  a  strong  dependence  between  the  source-drain  current  and  the  gate  electric  field.  These  results  illustrate  that 
substituents  and  symmetry  considerations  can  provide  considerable  control  of  the  behavior  of  molecular  transistors. 
Further  discussion  of  the  results  can  be  found  in  the  original  papers.12'14 

TEMPERATURE  EFFECTS  IN  MOLECULAR  DEVICES 

Chen  et  al.3  recently  reported  I-V  characteristics  of  molecules  consisting  of  chains  of  three  benzene  rings  with 
ligands  attached  at  various  places.  The  most  interesting  result  is  large  negative  differential  resistance  evinced  by  a 
relatively  sharp  spike  in  the  I-V  characteristic.  The  spike  is  found  to  broaden  and  shift  on  the  voltage  axis  with  increasing 
temperature.  The  shift,  by  about  1  V,  is  very  unusual.  In  semiconductor  nanostructures  resonant  peaks  have  been  found 
to  broaden  (by  standard  electron-phonon  interactions)  but  they  never  shift  appreciably. 

Calculations  for  three-benzene-ring  molecules  are  not  practical,  but  we  pursued  the  question  by  calculations  for  a 
single  benzene  ring  with  an  N02  ligand  replacing  one  of  the  hydrogen  atoms.13  We  found  that  the  energy  levels  of  the 
ligand  move  substantially  with  increasing  voltage  and  push  the  n  levels  into  the  active  window.  Thus  the  main  peak  in 
the  current  arises  primarily  from  tc  electrons  instead  of  n*  electrons.  We  then  explored  the  effect  of  rotating  the  ligand. 
We  found  that  a  rotation  by  90°  shifts  the  peak  to  lower  voltage  by  almost  1  V,  in  agreement  with  observations  (Fig.  4). 
The  interpretation  is  that  higher  temperatures  excite  the  rotational  modes  of  the  ligand.  Calculations  of  the  total  energy  of 
the  molecule  as  a  function  of  ligand  rotation  show  that  the  effective  rotational  quantum  of  energy  is  only  3  meV.  Thus, 
even  at  the  relatively  low  temperatures  of  the  experiment,  a  large  number  of  quanta  are  excited,  making  the  ligand 
effectively  a  classical  rotor  that  spends  most  of  its  time  at  the  extrema  of  the  amplitude.  Ligand  rotation  is  of  course  a 
unique  phenomenon  of  the  molecular  world,  explaining  why  large  voltage  shifts  are  not  observed  in  semiconductor 
nanostructures. 

CONCLUSIONS 

The  calculations  summarized  above  show  that  theory  has  now  advanced  to  the  point  where  quantitative  predictions 
can  be  made  about  transport  in  single  molecules.  Such  calculations  are  expected  to  play  a  major  role  in  the  evolution  of 
molecular  electronics,  the  way  that  simple  drift-diffusion  calculations  of  current  in  semiconductor  structures  have  played 


in  the  evolution  of  silicon-based  microelectronics. 
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Fig.  1  Schematics  of  benzene  ring  configurations  between 
electrodes  and  I-V  characteristics.  Top  row:  Experiment  (Ref. 
5).  Middle  row:  Theory  with  S  atoms  only  between  benzene 
ring  and  macroscopic  electrodes  (Ref.  7).  Lower  row:  Theory 
with  Au  atoms  inserted  between  the  S  atoms  and  the 
macroscopic  electrodes  (Ref.  7). 
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Fig.  2  Schematic  of  the  setup  showing  the  left  and  right  quasi 
Fermi  levels;  the  I-V  characteristic  as  in  Fig.  1,  and  the 
densities  of  states  discussed  in  the  text. 


Fig.  3  Schematic  of  three-terminal  devices  using  two 
different  molecules  as  discussed  in  the  text  and  the 
theoretical  I-V  characteristics  as  a  function  of  the  gate 
voltage  for  a  small  source-drain  voltage. 


Fig.  4  Experimental  data  by  Chen  et  al.  (Ref.  6)  at  two 
temperatures  and  theoretical  I-V  curves  for  two  different 
orientations  of  the  ligand.  See  text  for  more  details. 
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We  report  first-principles  calculations  of  the  current-voltage  (l-V)  characteristics  of  a  molecular  device 
and  compare  with  experiment.  We  find  that  the  shape  of  the  7-V  curve  is  largely  determined  by  the 
electronic  structure  of  the  molecule,  while  the  presence  of  single  atoms  at  the  molecule-electrode  interface 
play  a  key  role  in  determining  the  absolute  value  of  the  current  The  results  show  that  such  simulations 
would  be  useful  for  the  design  of  future  microelectronic  devices  for  which  the  Boltzmann-equation 
approach  is  no  longer  applicable. 

PACS  numbers:  73.40Jn,  73.40.Cg,  73.40.Gk,  85.65.+h 


Conventional  Si-based  microelectronics  is  likely  to 
reach  its  limit  of  miniaturization  in  the  next  10-15  years 
when  feature  lengths  shrink  below  100  nm.  The  main 
problem  is  the  onset  of  quantum  phenomena,  e.g.,  tunnel¬ 
ing,  that  would  make  scaled-down  conventional  devices 
inoperable.  Successor  technologies  currently  under  de¬ 
velopment,  such  as  tunneling  field-effect  transistors  and 
single-electron  transistors,  are  in  fact  based  on  quantum 
phenomena.  For  the  ultimate  miniaturization  below  nm, 
devices  made  from  single  molecules  are  currently  attract¬ 
ing  attention.  Prototypes  have  already  been  fabricated. 
Reed  et  al  [1]  reported  /-V  characteristics  of  single 
benzene-  1,4-dithio late  molecules.  Alivisatos  and  co¬ 
workers  [2]  reported  similar  7-V  characteristics  of  semi¬ 
conductor  and  metal  nanoclusters  between  gold  electrodes. 
Dekker  and  co-workers  [3]  reported  transistorlike  behavior 
in  carbon  nanotubes.  Similar  devices  have  been  demon¬ 
strated  by  Avouris  and  co-workers  using  single-walled 
and  multi  walled  carbon  nanotubes  [4]. 

Theoretical  modeling  played  a  key  role  in  the  invention 
of  the  transistor  and  in  the  subsequent  development  of  in¬ 
tegrated  circuits.  Device  modeling  continues  to  provide 
indispensable  input  to  circuit  modeling  for  designing  logic 
and  memory  chips  and  microprocessors.  It  is  based  on 
a  “semiclassical  approximation”  that  treats  electrons  and 
holes  as  classical  particles,  except  that  their  kinetic  ener¬ 
gies  are  determined  by  the  semiconductor  energy  bands, 
most  commonly  in  the  effective-mass  approximation.  In 
this  scheme,  transport  is  governed  by  Boltzmann’s  equa¬ 
tion.  Most  industrial  modeling  is  done  in  the  drift-diffusion 
approximation,  to  a  lesser  extent  in  the  higher-order  hydro- 
dynamic  approximation,  and,  at  times,  by  solving  Boltz¬ 
mann’s  equation  directly  by  Monte  Carlo  techniques  [5]. 
For  nanodevices,  however,  when  quantum  phenomena  are 
dominant,  the  semiclassical  Boltzmann ’s  equation  does  not 
apply  Quantum  mechanical  simulations  are  needed.  So 
far,  only  semiempirical  approaches  have  been  employed  to 
investigate  transport  in  molecular  systems,  providing  use¬ 
ful  insights  into  the  fundamental  mechanisms  [6-9]. 


In  this  Letter  we  report  first-principles  calculations  of 
the  7-V  characteristics  of  a  molecule.  In  this  method,  the 
electrostatic  potentials  through  the  molecule  and  at  the 
contacts  are  calculated  self-consistently  without  empiri¬ 
cal  adjustments.  We  report  calculations  for  the  benzene- 

1. 4- dithiolate  molecule  for  which  experimental  data  are 
available  [1].  We  find  that,  when  the  molecule  is  placed 
between  two  electrodes  made  of  an  ideal  metal  (homoge¬ 
neous  electron  gas  or  jellium  model  [10]),  the  shape  of  the 
7-V  characteristic  is  determined  by  the  electronic  structure 
of  the  molecule  in  contact  with  the  electrodes  and  in  the 
presence  of  the  external  electric  field.  The  shape  is  es¬ 
sentially  the  same  as  that  of  the  experimental  curve.  The 
absolute  magnitude  of  the  current,  however,  is  more  than 
2  orders  of  magnitude  larger  than  the  experimental  values. 
We  investigated  the  origins  of  this  discrepancy  and  found 
that  insertion  of  a  single  gold  atom  at  each  metal-molecule 
contact,  as  suggested  by  the  experimental  setup,  leaves  the 
shape  of  the  7-V  characteristic  globally  unchanged,  but  re¬ 
duces  the  absolute  value  of  the  current  by  more  than  an 
order  of  magnitude,  reflecting  reduced  coupling  between 
the  £  states  of  the  gold  and  the  p  orbitals  of  the  molecule, 
in  the  directions  parallel  to  the  electrode  surfaces.  Replace¬ 
ment  of  the  single  gold  atoms  by  aluminum  atoms,  which 
have  pz  orbitals  in  the  relevant  energy  region,  raises  the 
value  of  the  current  by  about  1  order  of  magnitude.  Other 
factors  that  can  affect  the  absolute  value  of  the  current  are 
discussed  later  in  this  paper.  These  results  show  that  such 
calculations  can  provide  valuable  quantitative  information 
to  help  design  molecular  devices. 

We  begin  with  the  study  of  the  7-V  characteristic  of 
a  single  benzene- 1,4-dithiol  molecule  between  two  ideal 
metallic  contacts.  It  is  well  known  that  when  the  benzene- 

1. 4- dithiol  molecule  is  adsorbed  on  gold  surfaces  the  H  of 
the  thiol  terminations  desorbs  and  the  sulfur  atoms  at  each 
end  bond  strongly  to  the  Au(lll)  surfaces  [11].  The  re¬ 
maining  molecule  ( benzene- 1, 4-dithiolate)  is  then  simply 
the  one  represented  in  Fig.  1,  where  we  show  the  contour 
plot  of  the  electronic  density  in  the  benzene  ring  plane.  We 
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FIG.  1 .  Contour  plot  of  the  electron  density  of  the  molecule 
described  in  the  text.  The  dots  represent  the  positions  of  the 
atoms.  The  lines  represent  the  position  of  the  model  metal 
surfaces. 


assume  that  the  molecule  stands  perpendicular  to  the  metal 
surfaces.  The  molecule  has  it  bonding  and  ir*  antibonding 
orbitals  formed  by  the  carbon  and  sulfur  p  orbitals  perpen¬ 
dicular  to  the  ring  plane  and  a  bonds  due  to  the  in-plane 
orbitals  of  the  atoms. 

We  computed  the  /-V  characteristic  using  the  method 
developed  in  Ref.  [12].  The  electron  density  of  the  jel- 
lium  electrodes  is  taken  equal  to  the  value  for  metallic 
gold  {rs  ~  3).  The  electron  wave  functions  are  computed 
by  solving  the  Lippman- Schwinger  equation  iteratively  to 
self-consistency  in  steady  state.  Exchange  and  correlation 
are  included  in  the  density-functional  formalism  within  the 
local-density  approximation  [13].  All  atomic  positions  are 
kept  fixed  at  their  equilibrium  values  in  the  free  molecule. 
The  current  is  computed  from  the  wave  functions  of  the 
electrode-molecule  system.  The  differential  conductance 
is  then  calculated  as  the  derivative  of  the  current  with  re¬ 
spect  to  the  external  bias.  Small  variations  in  the  atomic 
positions  (of  the  order  of  0. 1  A)  change  the  current  by  less 
than  1%. 

The  calculated  I-V  characteristic  is  shown  in  the  bottom 
panel  of  Fig.  2.  The  experimental  curve  is  also  shown 
for  comparison  in  Fig.  2.  It  is  clear  that  the  shapes  of 
the  two  curves  are  similar,  but  the  absolute  magnitude  of 
the  current  and  conductance  is  quite  different.  We  will 
first  discuss  the  origins  of  the  shape  and  then  address  the 
question  of  absolute  values. 

We  focus  on  three  distinct  regions  in  the  calculated  con¬ 
ductance  curve:  the  initial  rise  (from  zero  bias  to  about 
1  V),  the  first  peak  at  2.4  V,  and  the  second  peak  at  4.4  V. 
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FIG.  2.  Top:  Experimental  7-V  characteristic  of  a  benzene- 
1, 4-dithiolate  molecule  measured  by  Reed  et  al  [1].  Bottom: 
Conductance  of  the  molecule  of  Fig.  1  as  a  function  of  the 
external  bias  applied  to  the  metallic  contacts. 


In  Fig.  3  we  show  the  calculated  density  of  states  of  the 
molecule  for  three  different  voltages,  namely,  0.01,  2.4, 
and  4.4  V  (the  density  of  states  shown  is  the  difference 
between  that  of  the  molecule -electrode  system  and  that  of 
the  electrodes  without  the  molecule).  The  zero  of  energy 
is  the  left  Fermi  level  so  that  the  right  Fermi  level  is  equal 


Energy  (eV) 

FIG.  3.  Difference  between  the  density  of  states  of  the  two 
semi-infinite  electrodes  with  and  without  the  benzene- 1,4- 
dithiolate  molecule  in  between,  for  three  different  voltages.  The 
left  Fermi  level  (£Fl)  has  been  chosen  as  the  zero  of  energy. 
The  labels  £fr  correspond  to  the  energy  position  of  the  right 
Fermi  levels.  The  three  curves  correspond  to  the  bias  voltages 
indicated. 
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to  the  external  bias  (in  Fig.  3  the  labels  correspond  to 
the  2,4  and  4.4  V  values  of  the  bias). 

The  initial  slow  rise  of  the  conductance,  which  is  also 
present  in  the  experimental  data,  represents  basically 
Ohmic  behavior.  Figure  3  (top  curve)  shows  that  the 
molecule  has  a  small  but  relatively  smooth  density  of 
states  through  which  current  can  flow.  The  cr  and  7 r 
bonding  states  lie  several  eV  below  the  Fermi  energy  so 
that  a  considerable  bias  is  needed  to  produce  nonlinear 
behavior. 

After  the  initial  increase  of  the  current  with  increasing 
bias,  a  first  conductance  peak  and  subsequent  valley  are 
observed.  This  is  also  observed  experimentally,  though  at 
somewhat  lower  voltages  and  with  a  smaller  peak-to-valley 
ratio.  The  peak  and  the  valley  are  due  to  resonant  tunnel¬ 
ing  through  7T*  antibonding  states,  which  are  the  first  to 
appear  in  the  window  of  energy  between  the  right  and  left 
Fermi  levels  (see  Fig.  3,  middle  curve).  Inclusion  of  the 
vibrational  coupling  would  smear  out  the  peak  and  gener¬ 
ate  a  lower  peak-to-valley  ratio  in  the  /-V  characteristic, 
bringing  closer  the  agreement  between  theory  and  experi¬ 
ment.  Note  also  that  the  bonding  a  and  it  states  are  altered 
significantly  by  the  bias  simply  because  the  different  atoms 
of  the  benzene  ring  are  at  different  potentials. 

Increasing  the  bias  further,  a  second  I-V  peak  is  found 
due  to  resonant  tunneling  with  7 r  bonding  states.  From 
Fig.  3  it  is  evident  that  the  7 t*  antibonding  states  are  still 
present  in  the  energy  region  between  the  right  and  left 
Fermi  levels  but  with  a  lower  peak  ratio  with  respect  to  the 
7 r  bonding  states.  The  discrepancy  between  the  theoretical 
and  experimental  peak  positions  is  consistent  with  known 
limitations  of  the  local-density  approximation  [14]. 

The  different  scattering  processes  that  occur  for  differ¬ 
ent  biases  can  be  studied  in  more  detail  by  looking  at  the 
local  density  of  states  along  the  direction  of  current  flow. 
The  local  density  of  states  integrated  between  the  left  and 
right  Fermi  levels  for  small  bias  is  plotted  in  Fig.  4(A). 
As  in  the  case  of  the  total  density  of  states,  we  plot  here 
the  difference  between  the  local  density  of  states  of  the 
molecule-electrode  system  and  that  of  the  electrodes 
without  the  molecule.  It  is  evident  that  the  sulfur-to-metal 
contact  has  a  very  low  density  of  states,  in  effect  constitut¬ 
ing  a  barrier  through  which  electrons  must  tunnel.  After 
tunneling  through  this  barrier,  the  electrons  encounter 
a  region  of  relatively  smooth,  higher  density  of  states 
corresponding  to  the  states  of  the  molecule.  For  biases 
in  the  neighborhood  of  2.5  V  [Fig.  4(B)],  the  electrons 
find  a  much  larger  density  of  states  originating  from  the 
7T*  antibonding  states  of  the  molecule  (resonant  tunneling 
condition).  The  antibonding  region  is  between  the  middle 
carbons  of  the  molecule.  Finally,  with  a  further  increase 
of  the  bias,  to  the  range  of  5  V,  resonant  tunneling  occurs 
through  bonding  states  of  the  molecule  [see  Fig.  4(C)]. 

We  now  turn  to  the  absolute  value  of  the  current  which 
differs  by  more  than  2  orders  of  magnitude  from  the  experi¬ 
mental  value.  We  attribute  part  of  this  discrepancy  to  the 


s  c  c  c  c  s 


FIG.  4.  Local  density  of  states  difference  between  that  of  the 
molecule-electrode  system  and  that  of  the  electrodes  without  the 
molecule,  integrated  between  left  and  right  Fermi  levels.  Volt¬ 
age  region  (A)  corresponds  to  the  linear  response  regime.  Re¬ 
gion  (B)  corresponds  to  the  first  resonant  tunneling.  Region  (C) 
corresponds  to  the  second  resonant  tunneling  bias  conditions. 

contact  geometry  and  chemistiy.  The  experiment  by  Reed 
et  al.  [1]  suggests  that  the  contacts  are  atomically  termi¬ 
nated,  which  means  that  the  sulfur  atoms  are  attached  to 
single  Au  atoms.  We  investigated  the  effect  of  such  a  pos¬ 
sibility  on  the  I-V  characteristic  by  introducing  a  single 
gold  atom  between  the  sulfur  and  the  model  metal  sur¬ 
face  at  each  contact.  The  shape  of  the  I-V  characteris¬ 
tic  (Fig.  5)  remains  essentially  the  same,  but  the  absolute 
value  of  the  conductance  decreases  by  almost  2  orders  of 
magnitude.  This  reduction  is  partly  due  to  the  constriction 
resistance  generated  by  the  geometry  of  the  contact  and 
partly  due  to  the  chemistry  of  the  contact  between  the  gold 
and  the  sulfur  atoms.  The  Au  atoms  contribute  s  states  at 
the  Fermi  level  while  the  sulfur  atoms  that  attach  to  them 
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FIG.  5.  Conductance  of  the  molecule  of  Fig.  1  with  one  Au 
atom  between  the  model  metal  surface  and  the  sulfur  for  each 
contact  as  a  function  of  the  external  bias  applied  to  the  metallic 
contacts. 


contribute  p  states.  Only  the  p  states  that  are  perpendicu¬ 
lar  to  the  electrode  surfaces  can  couple,  by  symmetry,  to 
the  s  states  of  the  gold  atoms,  forming  <r  bonds.  The  p 
states  of  the  sulfur  atoms  that  are  parallel  to  the  metal  sur¬ 
faces  do  not  couple  to  the  gold  ,9  states,  thus  breaking  the  tt 
scattering  channel  that  connects  the  electrode  to  the  main 
benzene  ring.  The  gain  in  resistance  is  thus  mainly  due 
to  the  type  of  contacts.  This  conclusion  is  confirmed  by 
a  calculation  where  the  single  gold  atoms  at  the  contacts 
are  replaced  by  aluminum  atoms.  For  a  bias  of  0.01  V  the 
resistance  changes  from  96  to  2.9  Mfl.  The  A1  atoms  can 
now  contribute  with  p  orbitals  at  the  Fermi  level  that  are 
parallel  to  the  electrode  surface  and  form  7 r  states  with  the 
p  orbitals  of  the  sulfur  atoms  similarly  oriented.  Finally, 
we  performed  a  test  calculation  by  positioning  the  S  atom 
in  front  of  the  center  of  a  triangular  pad  of  three  gold  atoms 
on  each  electrode  surface.  This  simulates  a  S  atom  on  top 
of  a  Au(l  11)  surface.  The  calculated  resistance  is  nearly 
the  same  as  the  one  for  the  sulfur  attached  to  the  model 
metal.  The  reason  for  this  is  that  the  s  states  of  each  Au 
atom  in  the  pad  form  a  state  of  s  symmetry  and  two  states 
of  p  symmetry  parallel  to  the  surface  of  the  electrodes  that 
can  couple  to  the  corresponding  sulfur  p  state. 

The  above  calculations  show  that  the  absolute  magni¬ 
tude  of  the  current  is  an  extremely  sensitive  function  of 
the  contact  geometry  and  chemistry.  For  the  case  at  hand, 
the  discrepancy  with  experiment  is  still  about  an  order  of 
magnitude.  There  are  of  course  additional  effects  that  can 
alter  the  value  of  the  current,  such  as  temperature  (hot  elec¬ 
trons  and  vibrational  coupling)  or  local  disorder  in  the  Au 
metal  near  the  contacts  that  can  produce  electron  localiza¬ 
tion  [15].  The  latter  effect  may  arise  from  the  breaking  of 
the  gold  wire  when  contacts  to  the  molecule  are  made  [1]. 
We  further  note  that  the  experimental  measurements  have 
an  uncertainty  of  at  least  a  factor  of  2  [1].  Additional  theo¬ 
retical  and  experimental  work  is  clearly  needed  to  further 
clarify  these  issues. 


In  conclusion,  we  have  simulated  from  first-principles 
the  I-V  characteristics  of  a  molecular  device  and  found 
that  the  atomic-scale  contact  geometry  and  chemistry  play 
a  critical  role  in  the  absolute  value  of  the  conductance. 
This  result  sheds  new  light  on  the  transport  properties  of 
molecular  devices  and  shows  that  such  calculations  would 
be  useful  in  designing  devices  and  engineering  contacts  for 
future  nanotechnology. 
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Experiments  and  theory  have  so  far  demonstrated  that  single  molecules  can  form  the  core  of  a 
two-terminal  device.  Here  we  report  first-principles  calculations  of  transport  through  a  benzene- 1, 

4-ditlholate  molecule  with  a  third  capacitive  terminal  (gate).  We  find  that  the  resistance  of  the 
molecule  rises  from  its  zero-gate-bias  value  to  a  value  roughly  equal  to  the  quantum  of  resistance 
(12.9  kfl)  when  resonant  tunneling  through  the  7 r*  antibonding  orbitals  occurs.  ©  2000  American 
Institute  of  Physics.  [S0003-695 1(00)04023-7] 


Using  molecules  as  the  active  components  of  devices 
was  proposed  a  number  of  years  ago. 1  In  recent  years,  these 
ideas  have  been  revised,  extended,  and  realized.  Molecular 
rectification  was  demonstrated  in  1993.2  More  recently,  Reed 
etal  investigated  the  benzene- 1,  4-dithiol  rings  as  possible 
switching  devices  in  molecular  electronics  and  found  their 
current- voltage  (/-  V)  characteristics  to  be  highly 
reproducible.*5  Alivisatos  and  co-workers4  reported  similar 
/-  V  characteristics  of  semiconductor  and  metal  nanoclusters 
between  gold  electrodes.  Dekker  and  co-workers5  reported 
transistor-like  behavior  in  carbon  nanotubes.  Similar  devices 
have  been  demonstrated  by  Avouris  and  co-workers  using 
single-  and  multiwalled  carbon  nanotubes.6  Apart  from  the 
few  reports  on  transistor  operation  in  nanotubes,  three- 
terminal  geometries  in  molecule-based  devices  have  not  been 
investigated  because  of  the  difficulty  of  realizing  a  gate  ter¬ 
minal  at  such  short  length  scales.  However,  a  three-terminal 
device  is  ultimately  the  desired  device  for  many  of  the  ap¬ 
plications  of  molecules  in  electronics.  At  a  fixed  small 
source-drain  bias,  the  gate  voltage  must  be  able  to  amplify 
the  current  by  orders  of  magnitude  in  order  for  the  device  to 
be  a  possible  alternative  to  conventional  metal-oxide- 
semiconductor  field-effect  transistors. 

On  the  theoretical  side,  the  transport  properties  of  small 
molecules  cannot  be  modeled  by  solving  Boltzmann’s  equa¬ 
tion  as  is  done  in  conventional  electronics.7  Transport  prop¬ 
erties  must  be  calculated  directly  from  electron  wave  func¬ 
tions  using  a  full  quantum-mechanical  treatment.  So  far.  only 
semiempirical  approaches  have  been  employed  to  investigate 
transport  in  molecular  systems,  providing  useful  insights  into 
the  fundamental  mechanisms.8" 10  We  recently  reported  the 
first  ab  initio  simulations  of  the  I—V  characteristics  of  a 
molecule. 1 1 

In  this  letter  we  report  parameter-free,  fully  quantum 
mechanical  calculations  of  a  three-terminal  molecular  de¬ 
vice,  namely  a  benzene-1,  4-dithiolate  molecule  attached  to 
two  electrodes  and  a  capacitive  gate.  The  device  is  formed 
by  replacing  two  opposite  H  atoms  of  a  benzene  ring  by  two 
sulfur  atoms.3,11  The  S  atoms  are  then  attached  to  two  gold 
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electrodes.  The  gate  is  introduced  as  a  capacitor  field  gener¬ 
ated  by  two  circular  charged  disks  at  a  certain  distance  from 
each  other.  The  axis  of  the  capacitor  is  perpendicular  to  the 
transport  direction.  The  disks  are  kept  at  a  certain  potential 
difference  with  the  Fermi  level  on  one  disk  equal  to  the 
source  Fermi  level  (see  Fig.  1).  The  axis  of  the  cylindrical 
capacitor  lies  in  the  plane  of  the  benzene  ring.12  We  will 
show  that  the  benzene- 1,  4-dithiolate  molecule  exhibits  am¬ 
plification  when  a  polarization  field  is  applied  perpendicular 
to  the  transport  direction.  In  particular  we  show  that:  (a)  the 
molecule  behaves  as  a  resonant-tunneling  transistor  and  (b) 
no  charging  effect  occurs  in  the  molecule  because  the  elec¬ 
trons  do  not  spend  enough  time  in  the  device  to  prevent 
additional  charge  from  entering  the  molecule. 

The  benzene- 1,  4-dithiolate  molecule  has  been  studied 
experimentally  by  Reed  and  co-workers'5  in  a  two-terminal 
geometry.  The  corresponding  /-  V  characteristics  have  been 
theoretically  investigated  by  the  present  authors.11  It  was 
found  that  the  shape  of  the  /-  V  characteristic  is  determined 
by  the  electronic  structure  of  the  molecule,  as  modified  by 
the  interaction  with  the  electrodes  and  the  presence  of  the 
external  bias,  and  is  essentially  the  same  as  that  of  the  ex¬ 
perimental  curve.11  On  the  other  hand,  the  presence  of  single 
atoms  at  the  molecule-electrode  interface  controls  the  abso¬ 
lute  value  of  the  current  by  orders  of  magnitude.11  Amplifi¬ 
cation  by  employing  a  third  terminal  has  not  been  addressed 
either  experimentally  or  theoretically. 


FIG.  I.  Scheme  of  the  three- terminal  geometry  used  in  the  present  study. 
The  molecule  is  sandwiched  between  source  and  drain  electrodes  along  the 
direction  of  electronic  transport.  The  gate  electrodes  are  placed  perpendicu¬ 
lar  to  the  molecule  plane.  The  Fermi  level  on  one  gate  disk  equals  the  source 
Fermi  level  while  the  other  electrode  is  at  a  higher  potential  V0 . 


0003-6951/2000/76(23)/3448/3/$  17.00  3448  ©  2000  American  Institute  of  Physics 

Downloaded  25  dan  2001  to  129.59.117.185.  Redistribution  subject  to  A1P  copyright,  see  http://ojps.aip.org/aplo/aplcpyrts.html. 


Appi.  Phys.  Lett.,  Voi.  76,  No.  23,  5  June  2000 


Di  Ventra,  Pantelides,  and  Lang  3449 


Ec(vo!ts/A)  FtG.  3.  Difference  between  the  density  of  states  of  the  two  semi-infinite 

electrodes  with  and  without  the  benzene- 1 ,4-dithioiate  molecule  in  between, 
FIG.  2.  Conductance  of  the  molecule  of  Fig,  1  as  a  function  of  the  external  for  three  different  gate  fields  EG  (in  units  of  VIA).  The  left  Fermi  level  has 

gate  field.  The  source- drain  bias  is  0.01  Y\  been  chosen  as  the  zero  of  energy.  The  right  Fermi  level  is  at  10  mV. 

For  the  three-terminal  geometry,  we  computed  the  I-  V  electrodes  and  on  the  molecule  plane)  component  of  the  con- 

characteristic  using  the  method  developed  in  Ref.  13  (see  tinuum  states  couples  to  the  7 r*  states.15  This  is  the  equiva- 

also  Ref.  11).  Since  we  are  interested  in  the  effects  of  polar-  lent  of  the  usual  quadratic  Stark  shift  effect  in  atoms  and 

ization  on  transport,  we  choose  for  convenience  the  geom-  molecules,16  The  effect  is  less  pronounced  for  the  tt  states 

etry  for  which  the  molecule  is  attached  to  two  bare  metals  because  they  are  much  more  separated  in  energy  from  the 

represented  by  a  homogeneous  electron  gas  or  jeliium  continuum  states  of  higher  energy  than  are  the  7r*  orbitals. 

model.14  We  assume  that  the  molecule  stands  perpendicular  The  antibonding  states  thus  shift  in  energy  and  eventually 

to  the  metal  surfaces  (see  Fig.  1 ).  The  electron  density  of  the  enter  into  resonance  with  the  states  between  the  right  and  left 

jeliium  electrodes  is  taken  approximately  equal  to  the  value  Fermi  levels,  separated  by  10  meV  (middle  curve  of  Fig.  3). 

for  metallic  gold  (rf=3).  The  electron  wave  functions  are  Increasing  the  bias  further  (bottom  curve  of  Fig.  3),  the 

computed  by  solving  the  Lippman-Schwinger  equation  it-  resonant-tunneling  condition  is  lost  and  a  valley  in  the  I—V 

eratively  to  self-consistency  in  steady  state.13  Tlie  current  is  characteristic  is  observed.  Finally,  as  the  gate  bias  is  in- 

computed  from  the  wave  functions  of  the  electrode-  creased  further,  the  current  starts  to  increase  approximately 

molecule  system  in  the  presence  of  the  bias  and  gate  voltage.  linearly  with  the  gate  bias.  This  can  be  readily  understood 

The  differential  conductance  is  then  calculated  as  the  deriva-  within  the  Kubo-Greenwood  formalism:17  increasing  the 

tive  of  the  current  with  respect  to  the  external  bias. 1 1  Since  in  gate  bias,  quasifree  electron  states  with  increasing  energy 

practical  realizations  of  this  device  the  gate  could  be  of  dif-  enter  the  window  of  energy  between  the  right  and  left  Fermi 

ferent  form  and  size  we  discuss  the  results  in  terms  of  ap-  levels  and  can  thus  contribute  to  transport.  Assuming  the 

plied  gate  field  along  the  axis  of  the  capacitor.  electronic  effective  masses  are  approximately  the  same  at 

The  calculated  /-  V  characteristic  as  a  function  of  the  different  energies,  the  conductivity  (and  thus  the  current 

gate  bias  is  shown  in  Fig.  2  for  a  source-drain  voltage  dif-  since  the  source-drain  bias  is  only  10  mV)  increases  with 

ference  of  10  mV.  After  a  region  of  constant  resistance  the  the  square  of  the  density  of  states  at  the  energy  of  interest.17 

current  increases  with  the  gate  field,  reaches  a  maximum  For  the  electron  gas  of  these  quasifree  states,  the  density  of 

value  at  1.1  V/A,  then  decreases  at  about  1,5  V/A,  to  in-  states  can  be  assumed  to  vary  as  the  square  root  of  the  single 

crease  further  afterward  with  approximately  a  linear  depen-  particle  energy,  and  thus  the  current  varies  almost  linearly 

dence  on  the  gate  bias.  The  various  features  of  the  I-  V  curve  with  the  external  gate  bias. 

can  be  understood  by  looking  at  the  density  of  states  for  The  peak-to- valley  ratio  of  the  present  system  is  1.4, 

different  gate  voltages.  We  recall  here  that  the  molecule  has  which  is  probably  so  small  that  it  would  be  washed  out  by 

7 r  bonding  and  7 r*  antibonding  states  formed  by  the  carbon  the  vibrational  coupling  with  the  modes  of  the  molecule, 

and  sulfur  p  orbitals  perpendicular  to  the  ring  plane  and  cr  This  has  already  been  argued  in  the  two-terminal  geometry 

bonds  due  to  the  s  orbitals  and  in-plane  p  orbitals  of  the  case:11  the  peak-to-valley  ratios  observed  experimentally  in 

atoms.  the  present  system  are  considerably  smaller  than  theoreti- 

The  initial  slow  rise  of  the  conductance  represents  basi-  cally  predicted/’11 
cally  ohmic  behavior,  it  is  also  observed  experimentally  for  The  value  of  the  gate  field  at  which  resonant  tunneling 

the  two-terminal  geometry.3,11  Figure  3  (top  curve)  shows  occurs  («1  V/A)  seems  slightly  high  for  a  molecule  with  a 

that  the  molecule  has  a  small  but  relatively  smooth  density  of  nominal  length  of  8  A.  Two  observations  are  however  in 

states  through  which  current  can  flow.  The  tt  bonding  states  order:  (i)  we  found  in  previous  work11  that  the  theoretical 

lie  several  eV  below  the  Fermi  levels,  while  the  7 r*  antibond-  peak  of  transmission  due  to  7 r*  antibonding  states  occurs  at 

mg  states  are  nearly  I  eV  above  the  Fermi  levels.  higher  external  bias  than  the  experimental  one.  This  discrep- 

After  the  initial  increase  of  the  current  with  increasing  ancy  is  due  in  part  to  known  limitations  of  the  local  density 

gate  bias,  a  first  conductance  peak  and  subsequent  valley  are  approximation.18  (ii)  In  a  practical  realization  of  the  device, 

observed.  The  peak  and  the  valley  are  due  to  resonant  tun-  the  capacitance  field  would  certainly  leak  into  the  source  and 

neling  through  tt*  antibonding  states.  The  gate  field  couples  drain  electrodes,  providing  a  pocket  of  electrons  with  higher 

7 r*  states  with  the  continuum  of  states  at  higher  energy.  kinetic  energy  to  tunnel,  thus  reducing  the  gate  bias  value  at 

Since  the  gate  field  is  parallel  to  the  electrode  surfaces,  only  which  resonance  occurs.  Both  remarks  suggest  that  resonant 

the  odd  (with  respect  to  an  axis  parallel  to  the  surface  of  the  tunneling  and  amplification  can  occur  at  lower  gate  field. 
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The  resistance  at  zero  gate  bias  is  about  360  left.  At  LI 
V/A  its  value  changes  to  12.7  kft.  The  peak  value  is  close  to 
the  single-channel  resistance  of  12.9  kft  (spin  included).  At 
resonance,  the  lowest  7r*  antibonding  state  is  exactly  in  the 
narrow  window  of  energy  defined  by  the  right  and  left  Fermi 
levels,  providing  a  single  channel  for  electrons  to  tunnel  We 
stress  that  this  value  could  not  be  obtained  without  a  self- 
consistent  approach:  the  transmission  coefficient  of  the  chan¬ 
nel  is  not  necessarily  unity,10  and  depends  on  how  the  elec¬ 
tronic  charge  redistributes  in  the  molecule  and  at  the 
contacts.  The  amplification  (from  the  zero  gate  field  to  the 
gate  field  at  which  resonant  tunneling  occurs)  is  more  than  1 
order  of  magnitude.  The  measured  resistances  of  this 
molecule3  are  in  fact  much  larger  than  the  theoretical 
values.11  We  have  shown  that  these  resistances  are  partly  due 
to  the  chemistry  and  geometry  of  the  contacts  and  do  not 
affect  considerably  the  electronic  structure  of  the  molecule- 
electrode  system.11  Much  larger  amplification  should,  there¬ 
fore,  be  obtained  in  a  practical  implementation  of  the  device, 
since  the  value  of  the  resistance  when  resonant  tunneling 
occurs  is  close  to  (even  if  not  exactly  the  same  as)  the  single- 
channel  resistance. 

We  conclude  this  letter  by  suggesting  that  no  sizeable 
charging  effect  should  be  observed  in  the  I-  V  characteristics 
for  the  present  system.  Indeed,  the  electrons  do  not  spend 
enough  time  in  the  device  to  charge  it.  This  time  can  be 
estimated  by  calculating  the  delay  time  an  electron  spends  in 
the  scattering  region  whatever  the  outcome  of  the  scattering 
event  (the  electron  either  transmitted  or  reflected),39  and  is 
related  to  the  difference  of  density  of  states  AD(£)  of  the 
two  semi-infinite  electrodes  with  and  without  the  benzene-1, 
4-dithiolate  molecule  in  between  as  r=fnri\D(E).  At  small 
gate  fields  (where  Coulomb  blockade  could  potentially  be 
observed)  the  delay  time  is  about  2  fs.  In  this  time,  only  a 
charge  of  0.04|e|  is  present  in  the  molecular  region  between 
the  source  and  drain  electrodes  and  the  energy  region  be¬ 
tween  the  right  and  left  Fermi  levels.20  Due  to  the  short  delay 
time  and  die  small  extra  charge  present  in  the  molecular 
region,  we  suggest  that  transport  occurs  by  tunneling  and 
that  no  Coulomb  blockade  is  present. 


This  work  was  supported  in  part  by  Grant  No.  MDA972- 
98-1-0007  from  DARPA,  and  by  the  William  A.  and 
Nancy  F.  McMinn  Endowment  at  Vanderbilt  University.  The 
calculations  have  been  performed  on  the  IBM  SP2  at  the 
NPAC1  in  San  Diego. 
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The  conventional  Hellmann-Feynman  theorem  for  the  definition  of  forces  on  nuclei  is  not  directly  applicable 
to  quantum  time-dependent  and  transport  problems.  We  present  a  rigorous  derivation  of  a  general  Hellmann- 
Feynman-like  theorem  that  applies  to  all  quantum  mechanical  systems  and  reduces  to  well-known  results  for 
ground-state  problems.  It  provides  a  rigorous  definition  of  forces  in  time-dependent  and  transport  problems. 

Explicit  forms  of  Pulay-like  forces  are  derived  and  the  conditions  for  them  to  be  zero  are  identified.  A  practical 
scheme  for  ab  initio  calculations  of  current-induced  forces  is  described  and  the  study  of  the  transfer  of  a  Si 
atom  between  two  electrodes  is  presented  as  an  example. 


I.  INTRODUCTION 

The  Hellmann-Feynman  (HF)  theorem1  has  been  a  key 
ingredient  of  the  quantum  mechanical  treatment  of  forces 
acting  on  nuclei  in  molecules  and  solids.  In  turn,  these  forces 
are  the  key  ingredient  of  ab  initio  calculations  of  atomic- 
scale  structure  and  dynamics  in  materials  physics,  chemistry, 
and  molecular  biology.2"8  In  such  calculations,  the  electron 
system  is  kept  in  its  instantaneous  ground  state,  for  which  the 
traditional  formulations  of  the  HF  theorem  apply.  For  ex¬ 
ample,  for  several  decades,  theoretical  investigations  of  mol¬ 
ecules  and  chemical  reactions  have  relied  on  potential  energy 
surfaces  computed  in  this  fashion.  In  the  last  two  decades,  a 
similar  approach  has  been  the  basis  for  calculations  in  solids, 
e.g.,  surface  reconstruction,  phase  transformations,  defect 
configurations,  and  defect  reactions.  In  more  recent  years, 
fully  dynamical  calculations  (e.g.,  the  Car-Parrinello 
method2,3)  generally  assume  the  electron  system  remains  in 
its  instantaneous  ground  state. 

The  conventional  formulations  of  the  HF  theorem,2"8 
however,  do  not  apply  to  two  important  classes  of  problems 
in  which  the  electron  system  is  not  in  the  ground  state.  Both 
these  classes  are  emerging  frontiers  for  ab  initio  quantum 
mechanical  calculations. 

(a)  Molecules  and  solids  in  time-dependent  external  fields, 
e.g.,  radiation.  In  such  systems,  the  external  field  generally 
induces  electronic  excitations  and  nuclear  motions.  The  ad¬ 
vent  of  powerful  free-electron  lasers  and  of  tabletop  lasers 
that  deliver  intense  ultrashort  pulses  has  renewed  interest  in 
the  photodissociation  of  molecules  and  other  ultrafast 
reactions,9  the  generation  of  high  harmonics,  including 
x  rays,  from  atoms  in  intense  infrared  pulses,10,11  the  photo- 
induced  desorption  of  atoms  and  molecules  from  solid 
surfaces,12  and  other  phenomena,  especially  in  the  nonlinear 
regime. 

(b)  Transport  in  nanostructures  and  molecules.  The  fabri¬ 
cation  of  electronic  devices  using  nanoparticles  or 
molecules13,14  and  the  use  of  a  scanning-tunneling  micro- 

PRB  61 


scope  to  create  atomic-scale  structures  on  crystal  surfaces15 
represent  transport  problems  for  which  the  traditional  Boltz¬ 
mann  equation  approach  does  not  apply.16  Fully  quantum 
mechanical  calculations  of  transport  properties  are  needed, 
including  calculations  of  current-induced  atomic  rearrange¬ 
ments.  The  macroscopic  manifestation  of  the  latter,  namely 
electromigration,  has  remained  an  open  question  through  the 
years17,18  and  is  ripe  for  direct  first-principles  calculations. 

For  both  classes  of  problems,  the  conventional  formula¬ 
tions  of  the  HF  theorem  and  its  generalizations  are  not  ap¬ 
plicable  and  the  proper  definition  of  the  force  on  nuclei  has 
remained  unsettled. 10,1 8-20  Some  attempts  to  extend  the  HF 
theorem  to  nonstationary  states  have  been  proposed22  but,  as 
we  show  later  in  the  paper,  they  rely  on  an  incorrect  physical 
assumption. 

In  this  paper,  we  will  first  give  a  concise  statement  of  the 
HF  theorem  for  eigenstate  and  variational  ground-state  prob¬ 
lems  as  a  point  of  reference  and  identify  the  precise  obstacles 
that  have  left  the  problem  unsettled  for  the  time-dependent 
and  transport  problems.  We  will  then  present  a  general  form 
of  the  HF  theorem  that  is  applicable  to  all  quantum  mechani¬ 
cal  systems  with  square-integrable  wave  functions.  The 
proper  definition  of  forces  and  their  implementation  in  prac¬ 
tical  calculations  will  then  become  clear.  We  will  show  that 
this  general  theorem  and  its  practical  implementations  reduce 
to  the  classic  results  for  systems  in  the  ground  state.  An 
example  will  be  provided  for  the  forces  acting  on  a  Si  atom 
between  two  electrodes  when  steady-state  current  flows. 


n.  HF  THEOREM  FOR  GROUND  STATE 

Throughout  this  paper  we  will  present  derivations  using  a 
many-body  Hamiltonian  H  and  many-body  square-integrable 
wave  functions  ^  for  a  system  of  nuclei  and  electrons  with¬ 
out  any  particular  approximations.  Analogous  derivations 
can  easily  be  formulated  using  specific  approaches  to  the 
many-body  problem  such  as  Hartree-Fock  or  density- 
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functional  theory  (DFT).  The  explicit  equations  for  DFT  are 
given  in  the  Appendix. 

Given  the  eigenvalue  problem 

H9=E9,  (1) 

where  H  depends  parametrically  on  a  quantity  k  and  9  is 
square-integrable,  it  is  straightforward  to  show  by  direct  dif¬ 
ferentiation  that 

dE  I  SH  \  /  ,  \l  d9\ 

+  l-H-EVj  /  (t\V).  (2) 

In  view  of  Eq.  (1),  Eq.  (2)  reduces  to 

dE  I  3 H  \  /  , 

3r(*  n;*)/  ('W-  0) 

This  last  equation  represents  the  HF  theorem  for  exact  eigen¬ 
states.  If  the  quantity  \  is  a  given  degree  of  freedom  of  the 
system,  then  the  quantity  -  dEldk  given  by  Eq.  (3)  can  be 
interpreted  as  the  generalized  force  associated  with  it.  In  par¬ 
ticular,  if  the  nuclei  (or  ions)  are  treated  as  classical  particles 
and  \  is  the  position  vector  of  a  nucleus,  then  -dEldk  is 
the  classical  force  on  that  nucleus. 

The  practical  importance  of  the  HF  theorem  is  that  energy 
derivatives  are  difficult  to  compute  numerically,  whereas  the 
expression  on  the  right  in  Eq.  (3),  the  negative  of  which  is 
known  as  the  HF  force,  can  be  computed  efficiently.  How¬ 
ever,  it  was  recognized  by  early  computational  work  using 
atom-centered  basis  functions  that  the  HF  force  gave  mani¬ 
festly  wrong  results.5  The  origin  of  the  problem  is  most  suc¬ 
cinctly  illustrated  by  noting  that,  for  a  variational  calculation 
of  the  ground- state  energy  E ,  when  9  is  expanded  in  terms 
of  a  finite  basis  set,  one  no  longer  has  Eq.  (1)  but  the  more 
restrictive 

(9\H-E\9)  =  0,  (4) 

where  E  and  9  now  are  the  variational  energy  and  wave 
function,  respectively.  Equation  (4)  is  now  a  matrix  equa¬ 
tion.  As  a  result,  the  second  term  in  Eq.  (2)  is  no  longer  zero. 
Two  approaches  have  been  pursued  in  the  literature: 

(a)  One  requires  that  the  basis  set  is  such  that  the  extra 
term  in  Eq.  (2)  is  identically  zero.  This  requirement,  known 
as  Hurley’s  condition,21  is  satisfied  if  the  basis  functions  do 
not  depend  on  k  (e.g.,  when  the  parameters  k  are  nuclear 
position  vectors;  such  sets  are  known  as  “floating  sets”  in 
quantum  chemistry;  in  solid  state  applications,  the  com¬ 
monly  used  plane  waves  are  such  a  set).  The  condition  is 
also  satisfied  if  the  derivatives  of  the  basis  functions  with 
respect  to  k  are  themselves  part  of  the  basis  set.4  In  such  a 
case,  the  original  HF  theorem  is  satisfied.  In  order  to  prove 
these  statements,  one  writes  in  matrix  notation 

y  =  CXi  (5) 

where  {x}  is  the  basis  set  and  C  the  expansion  coefficients. 
From  the  stationarity  principle,  SE/S9  —  0,  we  get  SEISC 
=  0,  so  that  only  the  derivatives  dX/dk  survive  in  39/dk  in 
Eq.  (2). 


(b)  One  works  with  the  full  Eq.  (2).  The  usual  interpreta¬ 
tion  of  this  equation  is  that  -dE/dk  represents  the  exact 
force,  the  negative  of  the  first  term  on  the  right  is  the  HF 
force,  and  the  negative  of  the  last  term  is  a  correction  term 
often  referred  to  as  the  Pulay  force.5  This  interpretation  was 
backed  by  the  fact  that,  as  we  already  mentioned,  attempts  to 
use  the  HF  force  alone  with  atom-centered  basis  orbitals 
gave  manifestly  wrong  results.5  The  Pulay  force  is  routinely 
included  in  calculations  using  atom- centered  basis  functions 
and  is  in  fact  needed  to  produce  realistic  results.  It  is  viewed 
as  a  correction  to  the  HF  force  needed  to  make  the  calcula¬ 
tion  accurate  to  second  order.5,6 

From  this  perspective,  it  is  clear  why  the  problem  of 
forces  in  time-dependent  and  transport  problems  remains  un¬ 
solved:  it  is  not  obvious  that,  in  these  cases,  one  can  mean¬ 
ingfully  define  a  force  in  terms  of  an  energy  derivative. 
Moreover,  it  was  recognized  that  there  is  no  HF  theorem  that 
would  connect  the  energy  derivative  to  the  practical  HF 
force.18  It  has  been  argued,  however,  that,  even  in  the  ab¬ 
sence  of  a  HF  theorem,  the  HF  force  [right-hand  side  of  Eq. 
(3)]  remains  the  correct  quantum  mechanical  force  for  time- 
dependent  and  transport  problems  because  of  the  Ehrenfest 
theorem  that  relates  this  force  to  the  time  derivative  of  the 
expectation  value  of  the  (generalized)  momentum 
operator. 10,18  This  assertion  is  appealing,  but  has  not  been 
tested  and  runs  contrary  to  the  established  fact  that,  for 
variational  ground-state  calculations  using  atom-centered  ba¬ 
sis  sets,  the  HF  force  alone  is  inadequate  and  the  corrections 
introduced  by  the  Pulay  terms  are  not  negligible. 

m.  TIME-DEPENDENT  PROBLEM 

We  have  now  laid  the  groundwork  to  present  a  rigorous 
formulation  of  forces  for  all  quantum  mechanical  problems. 
The  first  step  is  to  recognize  that  the  only  correct  and  most 
general  quantum  mechanical  definition  of  a  force  is  the  ex¬ 
pectation  value  of  the  time  derivative  of  the  momentum 
operator.23  For  the  most  general,  time-dependent  quantum 
mechanical  problem  defined  by 

HQ  =  iytQ  (6) 

(we  use  units  such  that  fi- 1  and  we  use  <E>  to  denote  time- 
dependent  wave  functions),  the  force  for  a  given  degree  of 
freedom  k  is  defined  by23 


If  k  is  the  position  vector  of  a  given  nucleus,  then  the  force 
in  Eq.  (7)  is  simply  the  time  derivative  of  the  expectation 
value  of  the  momentum  operator  of  that  nucleus. 

If  we  specialize  this  definition  to  time-independent  prob¬ 
lems  for  which 

(t>(t)  =  e-iEt99  (8) 

we  get  immediately 


(9) 


PRB  61 


HELLMANN-FEYNMAN  THEOREM  AND  THE  DEFINITION  . . . 


16  209 


Two  most  important  observations  are  in  order:  first,  the 
reduction  of  the  force  definition  to  the  energy-derivative 
form  is  only  possible  for  time-independent  problems ;  and 
second,  for  time-independent  problems,  this  reduction  is 
valid  for  both  exact  wave  functions  and  wave  functions  ex¬ 
pressed  in  terms  of  a  finite  basis  set.  These  observations 
validate  rigorously  the  long-standing  assumption  that,  for 
time-independent  eigenstate  problems,  the  force  can  be  de¬ 
fined  in  terms  of  the  derivative  of  the  energy.  They  also 
make  clear  the  fact  that,  for  time-dependent  problems,  the 
rigorous  definition  of  the  force  is  Eq.  (7). 

We  now  address  the  issue  of  connecting  the  rigorous 
force  definition  to  the  HF  force  in  the  general  case.  By  direct 
differentiation  of  the  definition  (7)  and  use  of  Eq.  (6),  it  is 
straightforward  to  show  that,  for  exact  wave  functions 


d 

lJt 


3> 


(10) 


This  is  the  usual  Ehrenfest  theorem  that  can  be  found  in 
textbooks.23  In  view  of  Eqs.  (7)  and  (9),  we  see  immediately 
that  for  time-independent  problems,  the  Ehrenfest  theorem 
reduces  to  Eq.  (3),  namely  the  HF  theorem.  In  fact,  the 
Ehrenfest  theorem  has  been  referred  to  as  the  time-dependent 
version  of  the  HF  theorem.19  One  observation  is  crucial, 
however.  The  Ehrenfest  theorem  [Eq.  (10)]  is  only  valid  for 
exact  wave  functions.  When  the  wave  functions  are  ex¬ 
panded  in  terms  of  a  finite  basis  set,  as  in  most  computa¬ 
tional  work,  one  must  derive  the  corresponding  finite-set 
Ehrenfest-like  theorem .  We  note  that,  for  finite  basis  sets,  we 
no  longer  have  Eq.  (6),  but  the  more  restrictive 


or  if  their  derivatives  are  also  members  of  the  set.  In  order  to 
prove  this  statement,  it  is  necessary  to  invoke  the  fact  that  in 
time-dependent  problems,  the  action  S  defined  by 


Jtn 


<D 


H-i 


dt 


<S>)dt 


(13) 


is  stationary  for  a  given  d>,  namely  SSI  S^  =  0,  so  that  when 
we  expand  <E >—Cx,  we  have  SSI  SC = 0.  For  plane  waves 
and  “floating  sets,”  the  Pulay-like  terms  are  again  zero  and 
the  HF  force  is  the  exact  force.  Otherwise,  the  Pulay-like 
terms  must  be  included.  Note  also  that,  in  the  most  general 
case,  stationarity  of  the  action  does  not  imply  a  minimum 
principle,  indicating  that  the  accuracy  of  the  results  is  very 
sensitive  to  the  choice  of  the  approximate  wave  functions 
employed.10 

We  can  now  compare  directly  the  above  results  with  prior 
literature.  As  we  already  remarked,  the  assertion  by  Gross 
et  al 10  regarding  the  use  of  the  HF  force  is  rigorously  correct 
in  the  limit  of  a  complete  basis  set  and  remains  valid  when 
the  generalized  Hurley  condition  is  satisfied,  as  discussed 
previously  (see  also  the  Appendix).  Otherwise,  Pulay-like 
forces  must  be  included.  It  is  also  clear  from  the  above 
analysis  why  the  Bala  et  al  results22  are  not  correct:  the 
correct  definition  of  the  force  in  quantum  mechanics  is  the 
time  derivative  of  the  expectation  value  of  the  momentum 
operator  and  not  the  X  derivative  of  the  expectation  value  of 
the  Hamiltonian  on  time-dependent  wave  functions. 


IV.  STEADY-STATE  TRANSPORT  PROBLEM 


<D 


H-i- 


dt 


(ID 


which  is  a  matrix  equation.  By  direct  differentiation  and 
some  algebraic  manipulation,  we  then  get  the  finite-set 
Ehrenfest  theorem : 


d 

lJt 
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dH 

^X 
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H-i 
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This  is  the  central  result  of  this  paper  and  represents  the  most 
general  form  of  an  HF-like  theorem  that  is  applicable  to  all 
systems  and  allows  a  correct  and  unambiguous  definition  of 
forces  for  practical  implementation.  Note  that  this  approxi¬ 
mate  form  of  the  Ehrenfest  theorem  reduces  to  the  exact 
form,  Eq.  (10),  when  the  basis  is  complete,  simply  because 
Eq.  (6)  and  its  complex  conjugate  are  then  exact.  When  we 
specialize  the  finite-set  theorem  of  Eq.  (12)  to  the  time- 
independent  case,  we  immediately  get  Eq.  (2).  Thus,  as  is  the 
case  for  time-independent  problems,  we  find  that  when  the 
appropriate  form  of  the  Ehrenfest  theorem  is  used,  the  HF 
force  is  not  necessarily  the  correct  definition  of  force.  The 
extra  terms  in  Eq.  (12)  are  the  analog  of  the  Pulay  forces  for 
the  time-dependent  problem.  The  Pulay  forces  can  again  be 
made  zero  by  an  appropriate  choice  of  a  basis  set,  as  in  the 
time-independent  Hurley’s  condition.  In  other  words,  the  Pu¬ 
lay  forces  are  zero  if  the  basis  functions  do  not  depend  on  X 


We  turn  now  to  the  transport  problem.  The  general  for¬ 
malism  above  is  actually  valid  for  all  problems  describable 
by  square-integrable  wave  functions.  The  only  question, 
therefore,  is  one  of  normalization  of  wave  functions.  In  the 
eigenstate  and  variational  ground-state  problems  one  deals 
with  either  finite  systems  (molecules,  clusters)  or  with  infi¬ 
nite  solids  for  which  square-integrability  is  attained  through 
periodic  Bom- von  Karman  boundary  conditions.  For  trans¬ 
port  problems,  square-integrability  is  ensured  if  one  de¬ 
scribes  the  complete  circuit,  including,  e.g.,  the  battery  or 
generator.  In  the  most  general  case  of  time-dependent  power 
source,  the  above  time-dependent  formalism  is  completely 
valid.  For  practical  calculations,  however,  one  normally 
treats  a  system,  say  a  device  D,  in  contact  with  electrodes 
that  act  as  reservoirs  of  electrons  and  serve  merely  as  bound¬ 
ary  conditions  for  the  wave  function  of  D.  Square  integrabil- 
ity  can  then  be  assured  by  constructing  new  many-electron 
wave  functions  from  wave  packets  centered  at  individual 
one-electron  energies. 

For  time-independent  direct-current  (dc)  steady-state 
transport,  the  most  general  form  of  the  many-body  wave 
function  of  the  system  D  is  again  as  in  Eq.  (8),  where  now  E 
is  a  phase  factor  with  units  of  energy.  Using  this  form  in  Eq. 
(13)  we  find  that  the  action  becomes 

S«<®|ff-£|*)(r-f0).  (14) 

The  stationarity  of  the  action10  now  yields 

<5(<D|tf-£|d>)  =  0. 


(15) 
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The  equation  above  shows  that  the  dc  transport  problem  can 
be  rigorously  mapped  onto  a  variational  problem  where 
^|/7|Tjr)/^|ijf)  is  steady-state  energy  of  the  system. 

It  differs  from  the  ground-state  problem  only  in  the  form 
of  boundary  conditions  (closed  versus  open).  Again,  for 
plane  waves  and  “floating  sets/’  the  HF  force  represents  the 
total  force  whereas  for  other  sets  Pulay-like  forces  must  be 
included  as  in  Eq.  (2). 

So  far  our  derivations  were  carried  out  for  the  many-body 
Hamiltonian  H  and  wave  functions  W  or  4>.  For  practical 
implementations,  one  normally  separates  the  nuclear  and 
electronic  degrees  of  freedom  and  treats  the  nuclei  as  classi¬ 
cal  particles  (see  the  Appendix).  The  Hamiltonian  then  de¬ 
pends  parametrically  on  the  nuclear  position  vectors,  which 
can  be  treated  as  X’s  in  the  present  theory.  The  correspond¬ 
ing  forces  are  the  classical  forces  on  the  nuclei.  For  ground- 
state  and  time-dependent  problems  the  general  theory  above 
can  then  easily  be  cast  within  the  respective  Hartree-Fock  or 
density-functional  formulations  of  the  many-electron  prob¬ 
lem.  For  dc  transport,  Eq.  (15)  allows  us  to  conclude  that 
ground-state  density-functional  theory  is  again  applicable 
with  the  steady-state  energy  of  the  system  D  playing  the  role 
of  the  ground-state  energy.  This  result  provides  formal  justi¬ 
fication  for  the  use  of  one-electron  theories  for  transport  cal¬ 
culations.  We  note,  however,  that  for  the  transport  problem, 
a  practical  implementation  is  still  lacking  for  realistic  sys¬ 
tems.  We  conclude  the  paper  by  summarizing  the  main  ele¬ 
ments  of  such  a  practical  scheme. 

We  consider  the  system  D  and  two  (or  more)  electrodes 
(R  and  L  for  right  and  left).  Each  of  these  systems  can  be 
treated  separately  by  density  functional  theory.  Let  us  desig¬ 
nate  the  corresponding  single-particle  Hamiltonians  by  HD , 
Hr  ,  and  Hl  .  The  coupling  between  the  system  and  the  elec¬ 
trodes  can  then  be  included  by  defining  the  total  single¬ 
particle  Hamiltonian  as 

H=HD+HL+HR+rfLC+tfRC,  (16) 

where  Hllc  and  HRC  are  Hermitian  operators  to  be  deter¬ 
mined.  Only  the  Hamiltonians  H!LC,  HRC ,  and  HD  are  as¬ 
sumed  to  depend  on  the  position  of  the  atoms.  The  external 
electric  field  is  introduced  by  requiring  that,  far  from  the 
system,  the  Fermi-level  difference  between  the  left  and  right 
electrodes  is  equal  to  the  desired  value.20  As  we  showed 
above,  the  steady-state  energy  of  the  system  is  defined  in 
terms  of  the  ground-state  energy  functional  evaluated  with 
the  new  wave  functions  obeying  the  transport  boundary 
conditions.24 

The  most  effective  approach  to  a  self-consistent  determi¬ 
nation  of  H!lc,  Hrc,  and  *P  is  the  Green’s  function  method, 
where  HILC+HIRC+ffD  is  viewed  as  the  perturbation  to  the 
bare  electrodes.  So  far  this  problem  has  been  addressed  with¬ 
out  considering  the  problem  of  forces  by  several 
authors 20,25,26  A  practical  scheme  has  been  developed  by 
Lang20  by  assuming  that  the  electrodes  are  described  by 
“jellium,”  namely,  a  homogeneous  electron  gas  with  a 
smeared-out  compensating  positive  background.  Formula¬ 
tions  using  real  metals  are  not  well-suited  for  practical  cal¬ 
culations  because  of  ambiguities  in  ordering  the  wave  func¬ 
tions  by  energy.26  We  propose,  therefore,  the  following 
scheme  for  practical  calculations,  including  force  calcula¬ 


tions:  The  electrodes  can  be  approximated  by  jellium  only 
far  from  the  system-electrode  junctions  so  that  the  electrode 
wave  functions  can  be  easily  ordered  by  the  asymptotic  be¬ 
havior  as  plane  waves.  We  anticipate  that  only  a  few  layers 
of  metal  atoms  would  be  needed  attached  to  jellium  con- 
tinua.  For  the  actual  calculations,  the  metal  layers  are  viewed 
as  part  of  the  system  in  the  cavity  between  jellium  elec¬ 
trodes.  Lang’s  formalism  can  then  be  used  to  determine  H 
and  'P  self-consistently,  and  current-induced  forces  on  atoms 
can  be  computed  as  in  variational  problems  in  accordance 
with  the  theorem  proven  in  the  present  paper.27  The  spurious 
resistance  introduced  by  the  jellium-metal  interface  can  be 
calculated  separately  and  subtracted.  We  note  that  the  dis¬ 
crete  as  well  as  the  continuous  spectrum  of  the  whole  system 
can  be  included  in  the  formalism  above  in  a  straightforward 
way.  This  is  of  extreme  importance  in  such  phenomena  such 
as  electromigration,  where  the  contributions  coming  from  the 
continuous  and  discrete  parts  of  the  spectrum  are  of  the  same 
order  of  magnitude.18 

We  conclude  with  a  practical  example  applying  the  for¬ 
mulation  above  to  the  problem  of  a  Si  atom  between  two 
electrodes.28  The  two  electrodes  are  kept  at  a  distance  of  4.5 
A.  An  external  bias  is  applied  between  the  right  and  left 
electrode,  the  left  electrode  being  positive  with  respect  to  the 
right  electrode.  The  spectrum  of  the  present  system  has  a 
discrete  and  a  continuum  component.  For  the  single-particle 
wave  functions  in  the  discrete  part  of  the  spectrum  ipi9 
square-integrability  is  satisfied  because  the  fa's  are  local¬ 
ized.  For  each  energy  in  the  continuum  we  build  square- 
integrable  wave  functions  4/  in  an  energy  region  A,  as  we 
stated  above: 


where  A  is  a  normalization  constant  and  ^ s  are  single¬ 
particle  wave  functions  in  the  continuum.  As  in  Ref.  29  we 
choose  plane  waves  to  represent  the  Hilbert  space.30  Accord¬ 
ing  to  Eq.  (2),  the  Pulay-like  terms  are  thus  identically  zero. 
The  force  F  acting  on  the  atom  at  position  R  due  to  the 
electron  distribution  as  modified  by  the  external  bias  is  thus 


The  sum  and  integral  in  Eq.  (18)  include  spin  variable  too.  In 
the  present  case,  no  ion-ion  interaction  must  be  included. 
The  continuum  integration  a  covers  the  part  of  the  spectrum 
occupied  by  the  electrons  at  a  given  bias:  at  T=Q,  from  the 
bottom  of  the  conduction  band  on  the  left  to  the  quasi-Fermi- 
level  on  the  right.31 

In  the  case  at  hand,  the  force  is  directed  along  the  direc¬ 
tion  perpendicular  to  the  electrode  surfaces.  The  results  are 
plotted  in  Fig.  1  for  two  different  external  bias  conditions:  0 
V  and  3  V.28  At  zero  bias  the  Si  atom  has  two  stable  and  one 
metastable  configurations.  The  stable  configurations  corre¬ 
spond  to  the  atom  at  about  1.2  A  from  the  two  metals.  The 
metastable  configuration  corresponds  to  the  atom  between 
the  two  electrodes.  The  stable  configurations  at  zero  bias 
coincide  exactly  with  the  equilibrium  positions  obtainable 
from  standard  density-functional  total-energy  calculations.32 
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FIG.  1 .  Force  on  a  single  Si  atom  between  two  electrodes  sepa¬ 
rated  by  4.5  A  for  two  different  external  biases:  0  V  and  3  V. 
Positive  force  pushes  the  atom  to  the  right.  Zero  distance  corre¬ 
sponds  to  the  atom  between  the  two  electrodes. 

The  work  necessary  to  bring  the  atom  from  one  electrode  to 
the  other  can  be  calculated  from  the  curves  in  Fig.  1  as  the 
integral  between  the  metastable  and  stable  configurations.  At 
zero  bias  the  work  done  is  about  0.5  eV.  Applying  3  V  ac- 
cross  the  electrodes  (the  left  electrode  is  at  the  higher  bias), 
the  stable  positions  shift  slightly  towards  the  right  electrode 
due  to  a  transfer  of  charge  on  the  left  electrode.  Finally,  the 
activation  barrier  to  bring  the  atom  from  its  left  stable  con¬ 
figuration  to  the  new  metastable  position  decreases  to  about 
0.1  eV.  The  atom  is  essentially  free  to  jump  to  the  right 
electrode  when  the  bias  is  3  V. 

V.  CONCLUSIONS 


electron  system  with  charge  density  n(r)  in  an  external  po¬ 
tential  V(r)  can  be  written  as  a  functional  of  the  density33 


E[n]  =  Ts 


V(r)n(r)dr 


+ 


n(r)n(r') 

k-r'| 


drdrf  +  EXQ[n]. 


(Al) 


Ts  is  the  kinetic  energy  of  noninteracting  electrons,  and 
EM  is  the  exchange-correlation  energy.  The  quantum  me¬ 
chanical  problem  is  solved  with  variational  wave  functions. 
According  to  Eq.  (3),  the  electronic  contribution  to  the  force 
acting  on  a  nuclear  degree  of  freedom  X  is  thus 


SE( r)  dn{ r) 
Sn  S\ 


dr , 
(A2) 


where  SE(r)ISn  is  the  variation  of  the  functional  with  re¬ 
spect  to  the  density.  The  first  term  in  Eq.  (A2)  is  the  HF 
force,  while  the  second  is  the  Pulay-like  term.  If  the  ionic 
potentials  are  spherically  symmetric  and  their  non-Coulomb 
parts  do  not  overlap,  a  classical  electrostatic  term,  represent¬ 
ing  the  ion-ion  interaction,  must  be  added  to  Eq.  (A2)  to 
obtain  the  total  force  on  nuclei. 

In  time-dependent  DFT,  Eq.  (6)  decouples  into  a  set  of 
coupled  equations  for  all  electrons  ijjt  and  nuclei  <pi  (with 
charge  Zt)  under  external  time-dependent  potentials  V(r,  t) 
and  F(X,f),  respectively, 


We  presented  a  general  HF  theorem  that  applies  to  all 
quantum  mechanical  systems  and  allows  a  rigorous  defini¬ 
tion  of  forces  in  all  cases,  including  those  where  a  finite  basis 
set  is  used  to  represent  the  system  wave  functions.  In  the 
latter  case,  Pulay-like  forces  arise  that  can  be  set  to  zero  with 
a  particular  choice  of  the  basis  functions.  We  also  presented 
a  practical  scheme  for  transport  calculations  in  nanostruc¬ 
tures  that  includes  current-induced  forces  on  atoms.  The  lat¬ 
ter  scheme  is  particularly  important  nowadays  to  provide 
valuable  insight  into  current  effects  on  molecular  devices. 
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APPENDIX:  FORCES  IN  DFT 


where  7\  is  the  nuclei  kinetic  energy  and  £XcM(k>0  is  the 
nuclear  exchange-correlation  potential.  Note  that,  in  our  no¬ 
tation,  V(r ,t)  and  V(K,t)  contain  the  electron-ion  and  ion- 
ion  interactions,  respectively.  Assuming  the  nuclei  as  classi¬ 
cal  particles  with  nuclear  distribution  ^(X,*)- 
-  <5(X-X;(0)  and  applying  Eq.  (12)  to  the  nuclear  motion, 
the  force  on  the  nuclei  is 


™=-VX|(0K(A,0-Z, 


f  ”(rf)  ;  f  VM0”^ 


We  explicitly  write  in  this  appendix  the  expressions  of  the 
HF  forces  in  DFT.  We  recall  that  the  total  energy  of  a  many- 


The  last  term  of  Eq.  (A4)  is  the  Pulay-like  force  in  time- 
dependent  DFT. 
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Recent  experiments  found  an  unusual  temperature-induced  large  shift  in  the  resonant-tunneling  voltage 
of  certain  molecules.  We  report  first-principles  calculations  showing  that  such  behavior  can  be  caused 
by  the  excitation  of  rotational  modes  of  ligands.  These  modes  have  classical  characteristics,  i.e.,  the 
maximum  excursion  is  dominant,  while  at  the  same  time  they  have  a  significant  effect  on  the  energy 
levels  responsible  for  resonant  tunneling.  The  proposed  mechanism  of  ligand  rotations  is  unique  to 
molecules  and  accounts  for  the  fact  that  the  effect  is  not  seen  in  semiconductor  nanostructures. 
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Interest  in  the  transport  properties  of  single  molecules 
is  increasing  due  to  their  possible  use  as  electronic  com¬ 
ponents  [I].  New  synthetic  reactions  allow  the  realization 
of  new  and  diverse  molecules,  resulting  in  devices  with 
novel  and  unique  properties  [2].  For  example,  molecules 
have  been  shown  to  operate  as  Coulomb  blockade  struc¬ 
tures  [3],  diodes  [4],  or  switching  devices  with  high 
negative  differential  resistance  (NDR)  even  at  room  tem¬ 
perature  [5,6].  Theoretical  investigations  have  elucidated, 
e.g.,  the  role  of  contact  chemistry  and  geometry  on  the 
current- voltage  characteristics  of  single  molecules  [7], 
molecule-length  dependence  of  the  conductance  [8],  role 
of  the  injection  energy  at  the  contacts  [9],  three-terminal 
device  operation  [10],  and  stress-induced  modification  of 
transport  properties  [11]. 

Molecules  are  intrinsically  different  from  semiconduc¬ 
tor  nanostructures,  and  one  expects  that  some  of  their 
transport  properties  would  be  unique.  The  recent  report  by 
Chen  et  al.  [6]  confirms  that  this  is  the  case  with  the  tem¬ 
perature  dependence  of  the  resonant-tunneling  peak.  The 
molecules  studied  by  Chen  et  al.  are  chains  of  three  ben¬ 
zene  rings  with  one  or  more  ligands  attached  to  the  carbon 
atoms  of  the  central  ring.  An  unusually  sharp  peak  with 
an  enormous  NDR  has  been  observed  at  very  low  tem¬ 
peratures  [6].  At  higher  temperatures,  the  peak  broadens, 
as  expected,  but  it  also  shifts  to  substantially  lower  voltage 
[6].  The  large  shift  (~1  V)  is  highly  unusual  and  contrasts 
sharply  with  the  well-known  behavior  of  resonant  tunnel¬ 
ing  in  semiconductor  hetero structures  where  practically  no 
voltage  shifts  are  observed  [12,13].  Seminario  et  al  [14] 
have  proposed  that  the  voltage  shift  is  caused  by  a  rotation 
of  the  middle  carbon  ring  plus  charging  of  the  molecule 
by  single  extra  electrons.  The  suggestion  was  motivated  by 
ground-state  calculations  of  the  electronic  structure  of  the 
same  molecule  and  assumption  that  the  molecule’s  lowest- 
unoccupied  molecular  orbital  (LUMO — in  this  case  re¬ 
ferred  to  as  7 r*  orbital)  mediates  resonant  tunneling. 
Our  recent  first-principles  calculations  [7]  of  transport 
in  a  single  benzene  ring  have  shown,  however,  that  the 
electronic  structure  of  molecules  changes  substantially 
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as  a  function  of  applied  voltage ,  making  it  difficult  to 
infer  transport  properties  from  ground-state  electronic 
properties. 

In  this  Letter  we  report  first-principles  calculations  that 
explore  the  role  of  ligands  in  determining  the  transport 
properties  of  the  benzene-ring  family  of  molecules.  The 
study  was  motivated  by  the  fact  that  a  three-benzene-ring 
molecule  without  ligands  has  effectively  zero  conductivity 
for  all  voltages  (up  to  about  5  V)  and  insertion  of  ligands 
produces  significant  differences  in  the  shape  of  the  /-V 
curve  [6].  As  an  additional  indication  we  performed 
ground-state  density -functional  calculations  [15]  on 
the  three-ring  molecule  which  revealed  that  rotation  of 
the  middle  ring  produces  no  perceptible  change  in  the 
electronic  structure  of  the  molecule  in  the  vicinity  of 
the  highest-occupied  molecular  orbital  (HOMO)-LUMO 
states,  whereas  rotation  of  the  ligand  does.  Our  full  trans¬ 
port  calculations  for  a  single  benzene  ring  with  an  NO2 
ligand  find  essentially  the  same  behavior  as  that  exhibited 
by  the  three-ring  molecule  measured  by  Chen  et  al.  [6,16]: 
rotation  of  the  ligand — activated  by  temperature — causes 
a  substantial  shift  in  the  resonant-tunneling  voltage. 
This  large  voltage  shift  is  due  in  part  to  the  rotation 
properties  of  the  NO2  group  and  in  part  to  the  different 
symmetry  of  the  states  localized  on  the  NO2  group  with 
respect  to  the  rr  orbitals  of  the  carbon  ring  responsible 
for  transport.  In  particular,  unlike  the  7 t  orbitals  of 
the  carbon  ring,  the  orbitals  on  the  N02  group  have  a 
nonzero  dipole  moment:  they  have  a  lower  symmetry 
with  respect  to  reflection  through  a  plane  bisecting  the 
molecular  axis  so  that  they  are  shifted  in  energy  by  a 
field  applied  along  the  axis.  The  perturbation  of  the  states 
induced  by  the  external  electric  field  is  thus  first  order. 
At  increasing  temperature,  these  states  thus  “push”  the  rr 
orbital  to  reach  the  resonant  tunneling  condition  at  lower 
voltages.  The  results  show  that,  unlike  in  mesoscopic 
devices,  temperature  can  produce  different  effects  on 
resonant  tunneling  in  molecular  structures.  This  must 
be  taken  into  account  in  the  design  of  future  molecular 
devices. 
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FIG.  1.  Scheme  of  the  molecular  structure  investigated.  The 
structure  is  the  benzene- 1,4-dithiolate  molecule  where  a  H  atom 
is  replaced  by  a  NO2  group.  All  atoms  lie  on  the  plane  defined 
by  the  carbon  ring.  The  sulfurs  attach  to  ideal  metallic  leads. 

The  molecule  we  investigated  is  depicted  in  Fig.  1.  The 
two  sulfur  atoms  make  contact  to  gold  surfaces  that  we 
model  with  ideal  metals  (uniform  positive  background 
model  or  jellium  model  [17]).  The  interior  electron 
density  of  the  jellium  electrodes  is  taken  equal  to  the 
value  for  metallic  gold  ( rs  «  3).  We  computed  the  /-V 
characteristic  using  the  method  developed  in  Ref.  [18]. 
The  electron  wave  functions  are  computed  by  solving 
the  Lippman-Sch winger  equation  iteratively  to  self- 
consistency  in  steady  state.  Exchange  and  correlation  are 
included  in  the  density -functional  formalism  within  the 
local-density  approximation  [18].  The  current  is  computed 
from  the  wave  functions  of  the  electrode-molecule  system. 
At  zero  temperature,  the  atoms  of  the  NO2  group  lie  on 
the  same  plane  defined  by  the  six-carbon  ring  (see  Fig.  1). 

We  first  discuss  the  /-V  characteristics  at  zero  tempera¬ 
ture  and  then  the  effect  of  the  rotation  of  the  NO2  group. 
The  calculated  I-V  characteristic  of  the  molecule  is  plotted 
in  Fig.  2.  The  current  is  effectively  zero  for  voltages  up 
to  about  0.5  V.  The  current  then  increases  almost  linearly 
with  external  bias  until  a  peak  and  a  valley  are  visible  at 
about  3.8  and  4.2  V,  respectively.  This  behavior  can  be  un¬ 
derstood  by  looking  at  the  density  of  states  (DOS)  for  dif¬ 
ferent  external  voltages.  In  Fig.  3  we  show  the  calculated 
DOS  of  the  molecule  for  three  different  voltages,  namely, 
0.01,  3.8,  and  4.2  V  (the  DOS  shown  is  the  difference  be¬ 
tween  that  of  the  molecule-electrodes  system  and  that  of 
the  electrodes  without  the  molecule).  The  zero  of  energy 
is  the  left  Fermi  level  so  that  the  right  Fermi  level  is  equal 
to  the  external  bias  (in  Fig.  3  the  labels  correspond 
to  the  3.8  and  4.2  V  values  of  the  bias).  Contrary  to  the 
case  of  the  benzene- 1,4-dithiolate  molecule  (i.e.,  the  same 
molecule  with  the  N02  group  replaced  by  a  H)  the  7 t  bond¬ 
ing  orbital  of  the  carbon  ring  lies  close  to  the  left  Fermi 
level  (at  about  0.5  eV  below  the  left  Fermi  level  for  a  bias 
of  0.01  V — see  Fig.  3)  [7].  This  is  due  to  the  presence 


FIG.  2.  Theoretical  I-V  curve  of  the  molecular  structure  of 
Fig.  1  for  0°  and  90°  of  rotation  of  the  N02  group  with  respect 
to  the  carbon-ring  plane. 

of  the  hybridized  p  orbitals  of  nitrogen  and  oxygens  that 
push  the  7 t  orbital  higher  in  energy.  These  states  lie  about 
1.7  eV  below  the  7 t  orbital  for  a  0.01  V  bias  (indicated  as 
NO p  in  Fig.  3).  Therefore,  these  states  do  not  contribute 
directly  to  transport.  Increasing  the  bias,  the  7 r  orbital  en¬ 
ters  in  resonance  with  the  left  Fermi  level  and  a  peak  in 
the  current  occurs  at  this  resonant  tunneling  condition.  In¬ 
creasing  the  bias  further,  the  resonant  tunneling  condition 
is  no  longer  satisfied  and  a  valley  in  the  I-V  curve  is  ob¬ 
served  [19].  This  is  completely  different  from  the  case  of 
the  molecule  without  the  NO2  group  where  the  first  reso¬ 
nant  tunneling  condition  is  due  to  7 t*  antibonding  states 
[7].  We  note  that  the  resonant  tunneling  condition  at  zero 
temperature  occurs  at  3.8  V,  i.e.,  at  a  higher  voltage  than 
in  the  experiment  [6].  We  note,  however,  that  in  the  ex¬ 
periment  three  carbon  rings  form  the  entire  molecule.  In 
this  case,  the  rr  orbitals  of  the  outer  rings  will  push  the  7 r 
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FIG.  3.  Difference  between  the  density  of  states  of  the  two 
semi-infinite  electrodes  with  and  without  the  molecule  of  Fig.  1 
in  between,  for  three  different  voltages.  The  left  Fermi  level 
(Efl)  has  been  chosen  as  the  zero  of  energy.  The  labels  EFR 
correspond  to  the  energy  position  of  the  right  Fermi  levels.  The 
three  curves  correspond  to  the  bias  voltages  indicated. 
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FIG.  4.  Total  energy  of  the  isolated  molecule  of  Fig.  1  for 
different  angles  of  rotation  of  the  NO2  group  with  respect  to 
the  carbon-ring  plane. 

orbital  of  the  central  ring  even  closer  to  the  left  Fermi  level 
thus  allowing  resonant  tunneling  at  lower  voltage.  Over¬ 
all  the  molecule  behaves  as  a  molecular  resonant  tunneling 
device. 

Let  us  now  discuss  the  effects  of  temperature  on  the 
7-V  characteristics.  Increasing  the  temperature,  vibrational 
modes  can  be  excited.  However,  usual  optical-phonon 
scattering  cannot  account  for  the  large  shift  of  the  current 
peak  observed  experimentally  [13).  To  confirm  this,  we 
performed  a  set  of  calculations  of  the  current  at  external 
voltages  close  to  the  resonant  tunneling  condition  where 
a  “frozen  phonon”  has  been  allowed  in  the  structure.  We 
considered  the  highest  frequency  C-H  vibrational  mode, 
and  we  found  that  displacements  of  the  H  atoms  along  the 
C-H  bond  reduce  the  current  but  do  not  shift  the  voltage 
at  which  resonant  tunneling  occurs.  On  the  other  hand, 
the  rotation  of  the  NO2  group  with  respect  to  the  plane 
of  the  carbon  ring  produces  a  different  effect.  This  rota¬ 
tional  mode  can  easily  be  activated  by  temperature.  This 
is  demonstrated  in  Fig.  4  where  we  plot  the  total  energy 
of  the  free  molecule  for  different  rotation  angles  of  the 
NO2  group  with  respect  to  the  plane  of  the  carbon  ring 
[15].  The  total  energy  difference  between  0°  and  90°  of 
rotation  is  about  0.3  eV.  An  estimate  of  the  energy  level 
separation  for  this  vibrational  mode  can  be  obtained  by 
fitting  the  curve  of  Fig.  4  to  a  simple  harmonic  potential 
up  to  about  50°  of  rotation  [20].  The  energy  levels  are 
thus  En  —  ha)(n  4-  1/2),  n  =  0, 1, . . . ,  and  co  = 

Here  k  is  obtained  from  a  fit  to  a  harmonic  potential  and 
7  =  2Mo(doo/2)2  is  the  moment  of  inertia,  where  Mo  is 
the  mass  of  the  oxygen  and  doo  is  the  distance  between  the 
two  oxygens  of  the  NO2  group.  Taking  into  account  all  pa¬ 
rameters,  the  energy  separation  is  about  3  meV.  The  NO2 
group  can  thus  rotate  easily  at  room  temperature,  practi¬ 
cally  as  a  classical  rigid  rotator,  and  the  system  will  spend 
most  of  the  time  at  the  highest  degree  of  rotation  possible  at 
a  given  temperature.  This  rotation  produces  the  same  effect 
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FIG.  5.  Density  of  states  for  three  different  voltages  for  the 
molecule  with  the  N02  group  rotated  by  90°  with  respect  to  the 
carbon-ring  plane.  The  same  conventions  apply  as  in  Fig.  3. 

on  the  J-V  characteristics  as  observed  experimentally.  This 
is  shown  in  Fig.  2  for  the  90°  rotation.  The  peak-to-valley 
ratio  is  reduced,  and  the  peak  occurs  at  about  1  V  less  than 
at  zero  temperature.  The  reason  for  this  can  be  rational¬ 
ized  by  looking  at  the  DOS  for  different  external  voltages. 
In  Fig.  5  we  plot  the  calculated  DOS  of  the  molecule  with 
the  N02  group  rotated  by  90°  for  three  different  voltages, 
namely,  0.01,  2.8,  and  3.6  V.  The  same  conventions  as  in 
Fig.  3  apply.  The  rotation  of  the  N02  group  pushes  the 
NO-related  p  orbitals  very  close  to  the  ring  it  orbital  (the 
energy  separation  is  now  less  than  0.5  eV).  Because  of 
symmetry,  the  dipole  moment  of  the  tt  orbital  is  zero  while 
the  dipole  moment  of  the  NO-related  p  orbitals  is  non¬ 
vanishing.  Therefore,  by  applying  an  external  electric  field 
the  NO-related  p  orbitals  will  be  subject  to  a  linear  Stark 
shift  while  the  ring  7 r  orbital  shows  a  quadratic  Stark  shift 
[21].  Therefore,  with  increasing  bias,  the  NO  p  orbitals 
get  closer  to  the  left  Fermi  level  faster  than  the  7 r  orbital. 
Since,  in  the  present  case,  the  two  orbitals  are  very  close  in 
energy,  the  77  orbital  of  the  benzene  ring  will  be  “forced” 
to  reach  the  resonant  tunneling  condition  at  a  lower  bias 
due  to  the  presence  of  the  NO  p  orbitals.  The  effect  is 
sizable  since  at  90°  the  voltage  shift  for  the  peak  in  the  7-V 
curve  is  about  1  V  as  observed  in  the  experiment  [6]. 

In  conclusion,  we  have  shown  that  thermally  excited  ro¬ 
tational  modes  of  single  molecules,  like  rotation  of  ligands, 
can  have  dramatic  effects  on  their  7-V  characteristics.  The 
result  has  no  counterpart  in  usual  mesoscopic  physics  and 
must  be  taken  into  account  in  the  design  of  future  molecu¬ 
lar  devices. 
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We  report  first-principles  calculations  of  current-induced  forces  in  molecular  wires  for  which  experi¬ 
ments  are  available.  We  investigate,  as  an  example,  the  effect  of  current-induced  forces  on  a  benzene 
molecule  connected  to  two  bulk  electrodes  via  sulfur  end  groups.  We  find  that  the  molecule  twists 
around  an  axis  perpendicular  to  its  plane  and  undergoes  a  “breathing”  oscillation  at  resonant  tunneling 
via  antibonding  states.  However,  current-induced  forces  do  not  substantially  affect  the  absolute  value  of 
the  current  for  biases  as  high  as  5  V,  suggesting  that  molecular  wires  can  operate  at  very  large  electric 
fields  without  current-induced  breakdown. 
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The  phenomenon  of  atom  motion  due  to  current  flow 
(electromigration)  has  been  extensively  studied  in  the  past 
both  from  the  fundamental  standpoint  and  for  its  impor¬ 
tance  in  microelectronics  [1-4].  Most  recently,  a  new 
electronics  is  emerging  that  envisions  the  use  of  single 
molecules  or  molecular  wires  as  fundamental  components 
in  electronic  devices  [5].  For  instance,  it  has  been  demon¬ 
strated  that  molecules  can  operate  as  Coulomb  blockade 
structures  [6],  transistors  [7,8],  diodes  [9],  or  switching 
devices  with  high  negative  differential  resistance  even  at 
room  temperature  [10-12].  Since  electromigration  has 
been  a  major  concern  in  conventional  microelectronics  due 
to  current-induced  device  breakdown,  the  question  arises 
as  to  whether  current-induced  forces  may  present  a  severe 
limitation  to  the  development  of  molecular  electronics. 

It  was  recognized  in  early  theoretical  work  [1-4]  that 
current-induced  forces  on  a  given  physical  system  depend 
strongly  on  the  microscopic  details  of  the  self-consistent 
electric  field  that  is  created  upon  scattering  of  the  electrons., 
across  the  region  of  interest.  Self-consistency  in  the  calcu¬ 
lation  of  the  local  electric  field  with  the  correct  scattering 
boundary  conditions  is  thus  essential  to  have  meaningful 
quantitative  results  on  current-induced  forces. 

In  this  Letter,  we  report  first-principles  calculations  that 
explore  the  role  of  current-induced  forces  on  molecular 
wires,  and  their  role  in  weakening  chemical  bonds  at  the 
contacts  and  in  the  wire.  Considering  as  an  example  cur¬ 
rent  flow  in  a  benzene  molecule  connected  to  two  bulk  elec¬ 
trodes  via  sulfur  end  groups  [10],  we  extract  general  trends 
on  the  stability  of  molecular  wires  under  current  flow. 
The  molecular  structure  investigated  represents  a  proto¬ 
type  molecular  device  showing  nonlinear  transport  prop¬ 
erties  [10].  It  has  been  investigated  both  experimentally 
[10]  and  theoretically  [13-16]  without,  however,  address¬ 
ing  the  issue  of  current-induced  forces.  Since  we  are  in¬ 
terested  in  the  role  and  magnitude  of  these  forces,  and  the 
actual  experimental  contact  geometry  and  structure  of  the 
molecular  device  are  not  known,  we  focus  on  the  molecular 


PACS  numbers:  73.40.Jn,  73.40.Cg,  73.40.Gk,  85.65.  +h 

structure  depicted  in  Fig.  1.  The  effect  of  contact  chem¬ 
istry  and  geometry  in  the  current  has  been  discussed  in 
Ref.  [13].  We  find  that,  under  current  flow,  the  molecule 
twists  around  an  axis  perpendicular  to  its  plane  and  un¬ 
dergoes  a  “breathing”  oscillation  at  resonant  tunneling  via 
antibonding  states  [17].  However,  current-induced  forces 
do  not  substantially  affect  the  absolute  value  of  the  current 
up  to  biases  as  high  as  5  V.  This  is  a  remarkable  result 
for  a  molecule  of  nominal  length  of  only  8  A.  At  exter¬ 
nal  voltages  larger  than  5  V,  the  contact  that  is  depleted  of 
electrons  during  current  flow  weakens  considerably  with 
a  consequent  dramatic  reduction  of  the  current.  This  sug¬ 
gests  that  molecular  wires  can  operate  at  very  large  biases 
without  current-induced  breakdown,  in  contrast  to  recent 
findings  in  atomic  gold  wires  that  have  been  found  to  break 
at  biases  of  1  to  2  V  [18,19]. 

We  computed  the  I-V  characteristics  of  the  molecular 
structure  by  using  the  method  discussed  in  Ref.  [20].  The 
two  sulfur  atoms  of  the  molecule  (see  Fig.  1)  make  contact 
to  gold  surfaces  that  we  model  with  ideal  metals  (jellium 
model)  [20].  The  interior  electron  density  of  the  elec¬ 
trodes  is  taken  equal  to  the  value  for  metallic  gold  (rs  ~  3). 
The  electron  wave  functions  are  computed  by  solving  the 
Lippman-Schwinger  equation  iteratively  to  self-consistency 
in  steady  state.  Exchange  and  correlation  are  included  in 


FIG.  1.  Scheme  of  the  molecular  structure  investigated.  The 
structure  is  the  benzene- 1 ,4-dithiolate  molecule.  All  atoms  lie 
on  the  plane  defined  by  the  carbon  ring.  The  sulfurs  attach  to 
ideal  metallic  leads. 
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the  density-functional  formalism  within  the  local-density 
approximation  [20].  The  current  is  computed  from  the 
wave  functions  |^)  of  the  electrode-molecule  system.  The 
force  F  acting  on  a  given  atom  at  position  R  due  to 
the  electron  distribution  as  modified  by  the  external  bias 
is  given  by  the  Hellmann-Feynman  type  of  theorem  devel¬ 
oped  in  Ref.  [21]: 


F  = 


dEifa 


dH 

aR 


The  sum  and  integral  in  Eq.  (1)  include  spin  variables  also. 
The  first  term  on  the  right-hand  side  of  Eq.  (1)  is  the  usual 
Hellmann-Feynman  contribution  to  the  force  due  to  local¬ 
ized  electronic  states  |  ^/).  The  second  term  is  the  contri¬ 
bution  to  the  force  due  to  the  continuum  of  states  [21].  It  is 
calculated  by  constructing,  for  each  energy  in  the  contin¬ 
uum,  square-integrable  wave  functions  \ifr&)  in  an  energy 
region  A 


FIG.  2.  Theoretical  I-V  curve  of  the  molecular  structure  of 
Fig.  1  with  and  without  the  effect  of  current-induced  forces. 
The  vertical  arrows  indicate  the  onset  of  resonant  tunneling  via 
antibonding  (at  —2.4  V)  and  bonding  (at  —4.4  V)  states. 


=  A  (2) 

where  JA  is  a  normalization  constant  and  the  ij/'s  are 
single-particle  wave  functions  in  the  continuum,  solutions 
of  the  Lippmann-Schwinger  equation  [21].  The  continuum 
integration  a  covers  the  part  of  the  spectrum  occupied  by 
the  electrons  at  a  given  bias  [22].  Finally,  the  total  force 
on  the  atom  includes  a  trivial  ion-ion  interaction.  Starting 
from  a  given  atomic  configuration  (e.g.,  the  atoms  at  the 
equilibrium  experimental  atomic  positions),  we  calculate 
the  forces  acting  on  each  atom.  We  then  move  the  atoms 
according  to  the  gradient  of  these  forces  until  the  force  on 
each  atom  is  zero. 

We  now  discuss  the  main  results  obtained  using  the 
above  theoretical  approach.  The  relaxed  configuration 
of  the  system  at  zero  bias  consists  of  C-C  bond  lengths 
of  1.40  A,  C-H  bond  lengths  of  1.09  A,  C-S  bonds  of 
1.70  A,  and  S-jellium  surface  bond  length  of  1.00  A.  The 
latter  is  in  agreement  with  the  equilibrium  distance  of 
sulfur  adsorption  on  jellium  surfaces  [23],  and  the  other 
bond  lengths  are  in  good  agreement  with  the  experimental 
bond  lengths  in  isolated  benzene  molecules  and  thiophenol 
molecules  [24]. 

The  I-V  characteristic  of  the  molecular  structure  with 
and  without  the  effect  of  current-induced  forces  is  reported 
in  Fig.  2.  The  first  peak  in  conductance  (indicated  as  a 
vertical  arrow  in  Fig.  2)  occurs  at  about  2.4  V  and  is  due 
to  resonant  tunneling  via  7r*  antibonding  states  (see  also 
Ref.  [13]).  The  second  peak  at  about  4.4  V  (also  indicated 
as  a  vertical  arrow  in  Fig.  2)  is  due  to  resonant  tunneling 
via  7r  bonding  states  [13].  The  electron  transport  via  anti¬ 
bonding  states  corresponds  to  a  depletion  of  charge  in 
the  central  C-C  bonds,  and  an  accumulation  of  charge  in 
the  nearest  C-C  bonds  (see  Fig.  3a).  The  total  charge 
depleted  with  respect  to  the  zero  bias  condition  is  about 
0.05^  per  bond.  This  charge  is  nearly  completely  recov¬ 


ered  in  the  middle  bonds  when  the  resonant  tunneling  con¬ 
dition  is  lost.  In  particular,  at  4.4  V — corresponding  to 
resonant-tunneling  via  tt  bonding  states — the  charge  is 


s  c  c  c  c  s 


s  c  c  c  c  s 


FIG.  3.  Local  density  of  states  difference  between  that  of  the 
molecule-electrodes  system  and  that  of  the  electrodes  without 
the  molecule,  integrated  between  left  and  right  Fermi  levels  for 
a  bias  of  2.4  V  (A)  and  4.4  V  (B).  Solid  lines  correspond  to  the 
unrelaxed  geometry,  dashed  lines  to  the  relaxed  one.  The  vertical 
dotted  lines  correspond  to  the  unrelaxed  atomic  positions. 
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more  uniformly  distributed  across  all  bonds  in  the  benzene 
ring  (see  Fig.  3b). 

It  is  immediately  evident  from  Fig.  2  that  current- 
induced  forces  do  not  substantially  alter  the  absolute 
value  of  the  current  for  external  voltages  as  high  as  5  V. 
Nonetheless,  the  molecular  structure  undergoes  some 
structural  transformations  (see  below). 

The  bonding  and  antibonding  nature  of  the  states  in¬ 
volved  in  the  electron  scattering  has  a  dramatic  effect  on 
the  dynamics  of  the  molecule  under  current  flow.  This  is 
illustrated  in  Fig.  4  for  different  external  biases.  The  struc¬ 
ture  at  0  V  has  been  found  to  be  unstable  under  current  flow 
(see  Fig.  4):  With  increasing  bias,  the  mirror  symmetry 
with  respect  to  a  plane  perpendicular  to  both  the  benzene 
ring  plane  and  the  surface  of  the  electrodes  can  be  easily 
broken,  leading  to  a  slight  rotation  of  the  central  carbon 
ring  with  respect  to  an  axis  perpendicular  to  its  plane  [17]. 
At  the  same  time,  the  S-metal  bond  on  the  right  electrode 
weakens  due  to  the  transfer  of  charge  from  the  right  to  the 
left  electrode  (the  left  electrode  is  at  a  positive  bias  with 
respect  to  the  right  electrode,  see  Fig.  3a).  At  about  2.4  V, 
i.e.,  when  resonant  tunneling  via  antibonding  states  occurs, 
some  charge  depletes  from  the  central  C-C  bonds,  leading 
to  a  weakening  of  these  bonds.  Consequently,  these  bonds 
slightly  expand  forcing  the  remaining  C-C  and  C-S  bonds 
to  expand.  These  bonds  expand  on  average  0.05  A  while 
the  C-H  bonds  are  not  affected.  This  can  be  rationalized 
by  knowing  that  the  77*  states  of  the  molecule  are  formed 
only  by  carbon  and  sulfur  p  orbitals  that  are  perpendicular 
to  the  ring  plane,  while  the  C-H  bonds  are  cr-like.  Upon 
relaxation,  some  extra  charge  is  depleted  from  the  central 
C-C  bonds  with  consequent  redistribution  in  the  nearby 
bonds  (see  Fig.  3a,  dashed  line). 

Increasing  the  bias  further,  the  resonant-tunneling  con¬ 
dition  is  lost  and  charge  is  almost  completely  recovered  in 
the  central  C-C  bonds  (see  Fig.  3).  The  C-C  bonds  then 
contract  back  to  almost  their  length  at  0  V.  In  summary, 
the  molecule  undergoes  a  “breathing”  oscillation  when  the 
bias  is  scanned  across  the  first  resonant-tunneling  condi¬ 
tion.  This  “breathing”  oscillation  as  a  function  of  bias  is 
typical  only  of  resonant-tunneling  via  antibonding  states. 
Indeed,  increasing  the  bias  further,  until  resonant-tunneling 
via  bonding  states  is  satisfied  (—4.4  V),  no  bond-length  os¬ 
cillations  are  observed  due  to  a  more  uniform  distribution 
of  charge  across  the  central  C  ring  (see  Fig.  3b).  On  the 
other  hand,  with  increasing  bias,  the  central  ring  continues 
to  twist  with  respect  to  an  axis  perpendicular  to  its  plane, 
while  the  S-metal  bond  on  the  right  electrode  weakens  (see 

FIG.  4.  Structural  transformations  of  the  molecule  of  Fig.  1 
for  four  external  biases.  At  2.4  V,  the  C-C  bonds  of  the 
molecule  slightly  expand,  while  at  2.8  V,  the  same  bonds 
contract,  as  indicated  by  the  arrows.  Upon  symmetry  breaking, 
a  counterclockwise  (or,  equivalently,  a  clockwise)  rotation  of 
the  central  carbon  ring  with  respect  to  an  axis  perpendicular  to 
its  plane  is  observed.  The  left  electrode  is  at  a  positive  bias 
with  respect  to  the  right  electrode. 
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Fig.  3b,  dashed  line).  Remarkably,  the  absolute  value  of 
the  current  depends  only  weakly  on  the  dynamical  changes 
of  the  molecule  for  voltages  up  to  about  5  V  (see  Fig.  2). 
At  these  high  voltages,  however,  the  charge  transfer  from 
the  right  to  the  left  electrode  strongly  weakens  the  S-metal 
bond  on  the  right  electrode.  This  bond  expands  by  more 
than  0.3  A  for  biases  larger  than  5  V  (at  6  V  the  bond  ex¬ 
pansion  is  0.5  A),  therefore  behaving  as  an  extra  barrier 
for  electrons  to  tunnel  across  the  molecular  structure.  The 
sign  of  the  current-induced  forces  on  the  S  atom  of  this 
bond  coincide  with  the  sign  of  the  electron  current  [25]. 
Because  of  thermal  and  current  fluctuations,  and  the  fact 
that  the  S-metal  bond  has  been  weakened,  the  S-metal  bond 
distance  can  oscillate  at  these  high  voltages,  giving  rise  to 
small  oscillations  in  the  conductance  [26].  Small  conduc¬ 
tance  oscillations  are  indeed  observed  for  such  high  volt¬ 
ages  in  the  present  system  [10].  It  is  also  worth  noticing 
that  complete  fracture  of  this  bond  can  occur  at  these  high 
voltages  if  temperature  effects  are  taken  into  account  [19]. 

In  conclusion,  we  have  shown,  using  first-principles  cal¬ 
culations,  that  current-induced  forces  in  molecular  devices 
can  induce  unusual  dynamical  changes  in  the  structure  of 
the  molecules.  However,  we  have  found  that  the  absolute 
value  of  the  current  is  quite  unaffected  up  to  external  volt¬ 
ages  as  high  as  5  V,  in  contrast  to  the  case  of  atomic  gold 
wires  that  break  at  smaller  biases:  The  strong  a  bonds 
of  the  carbon-based  molecular  structures  make  them  more 
resistant  to  current-induced  forces  than  atomic  gold  wires. 
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